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We present a general scheme based on nonlinear response theory to calculate the expansion of 
correlation functions such as the pair-correlation function or the exchange-correlation hole of an 
inhomogeneous many-particle system in terms of density derivatives of arbitrary order. We further 
derive a consistency condition that is necessary for the existence of the gradient expansion. This 
condition is used to carry out an infinite summation of terms involving response functions up to 
infinite order from which it follows that the coefficient functions of the gradient expansion can be 
expressed in terms the local density profile rather than the background density around which the 
expansion is carried out. We apply the method to the calculation of the gradient expansion of the 
one-particle density matrix to second order in the density gradients and recover in an alternative 
manner the result of Gross and Dreizler (Z. Phys. A 302, 103 (1981)) which was derived using the 
Kirzhnits method. The nonlinear response method is more general and avoids the turning point 
problem of the Kirzhnits expansion. We further give a description of the exchange hole in momentum 
space and confirm the wave vector analysis of Langreth and Perdew (Phys. Rev. B21, 5469 (1980)) 
for this case. This is used derive that the second order gradient expansion of the system averaged 
exchange hole satisfies the hole sum rule and to calculate the gradient coefficient of the exchange 
energy without the need to regularize divergent integrals. 

PACS numbers: 31.15.ee, 71.15.Mb, 31.10,+z 



I. INTRODUCTION 

Since the groundbreaking work of Hohcnberg and 
Kohn [l[ we know that the external potential of any inho- 
mogeneous quantum many-body system is a functional of 
its ground state density n. This implies that the many- 
body ground state and hence any ground state 
expectation value 



as 



0,a 



A[n] = (*[n]|i|*[n]) 



(1) 



for any operator A is a functional of the density. This 
idea has inspired an enormous amount of work in a re- 
search field that is now known as density-functional the- 
ory. Density-functional theory 0-0] is currently one of 
the main approaches used in electronic structure the- 
ory. Over the years many extensions beyond its origi- 
nal formulation have been developed and currently it is 
widely applied in solid state physics and quantum chem- 
istry, both for ground state and time-dependent systems 
0. Especially a large activity has gone into finding ex- 
plicit expressions for the ground state energy E[n] as a 
functional of the density for many-electron systems. The 
ground state energy is obtained by minimizing this func- 
tional on an appropriate set of electronic densities, which 
is of great importance in determining structures and ge- 
ometries of molecules and solids. By use of the Kohn- 
Sham method @ this minimization problem is reduced 
to solving effective one-particle equations. The construc- 
tion of the effective potential in these equations requires 
the knowledge of the so-called exchange-correlation en- 
ergy functional E xc [n] . This functional can be expressed 



E*M = 5 J dvdv | r _ r ,| 



(2) 



where p xc (r'|r) is the coupling constant integrated 
exchange-correlation hole. The exchange-correlation hole 
has a physical interpretation as the difference between 
the conditional and the unconditional probability (which 
is simply the density n(r')) to find a electron at r', given 
that we know that there is an electron at r. Equation ([2|) 
is an important relation since it expresses the exchange- 
correlation energy in terms of a quantity that has a lo- 
cal physical interpretation and can be studied by accu- 
rate wave function and many-body methods or modeled 
based on physical intuition. It has therefore played an 
important role in the development of approximate den- 
sity functionals fsl— Tl^3| - One of the simplest and already 
very successful approximations has been the local density 
approximation (LDA) in which the exchange-correlation 
hole is taken from the homogeneous electron gas and ap- 
plied locally [|| to systems with arbitrary density profiles. 
The accuracy of the LDA has been considerably improved 
by means of the so-called generalized gradient approxi- 
mations (GGAs) [l5|. A large class of commonly used 
GGAs is based on the so-called real-space cutoff of the 
straightforward gradient expansion for the exchange and 
correlation hole. The first such functional for exchange 
was proposed by Perdew Q . He noted that the gradient 
expansion of the exchange hole to second order in the 
density gradients violates the negativity condition and 
the hole sum rule by which the exchange hole integrates 
to a deficit of one electron. By enforcing these physical 
constraints a density functional was obtained that greatly 



2 



improved on the LDA for binding energies of molecules. 
This procedure was later simplified [lfj using the fact 
that the exchange correlation energy ([2]) can be written 

as 



N f 
N = y J dy 



(pxc(y)) 



where N is the number of electrons in the system and 
(p X c(y)) = Jj J dr?i(r)p xc (r + y|r), (3) 

is the so-called system-averaged exchange-correlation 
hole. This averaged hole still obeys the negativity con- 
dition as well as the sum rule and can therefore be sub- 
jected to the real-space cut-off procedure. One advan- 
tage of using the system averaged holes was to reduce 
the order of the derivatives in the gradient expansion. 
Furthermore it also allowed for the real-space cut-off pro- 
cedure to be applied to the correlation hole since there is 
a known gradient expansion for the system averaged cor- 
relation hole in the random phase approximation (RPA) 
but no known expansion for the correlation hole itself. 
This was the basis for a GGA for the correlation energy 
SEMI- These GGAs relied heavily on the pioneering 
work by Langreth, Perdew and Mehl [l9l - |2~lj ] who per- 
formed a wave-vector decomposition of the system aver- 
aged hole of Eq. and calculated the Fourier transform 



(Pxc(k)) = / dy(p xc {y))e 



-iky 



from which the exchange correlation energy can be cal- 
culated as 



N 



We are not aware of any work that goes beyond refer- 
ences [l9l - |2~l| and improves on the straightforward gra- 
dient expansion of the system averaged correlation hole. 
The knowledge on the gradient expansion of the correla- 
tion hole is very limited indeed. As mentioned before, in 
contrast to the work on the exchange hole 0, @, there 
are no known expressions for the gradient expansion of 
the non-system averaged correlation hole. 
Part of the problem is that there has been no clear deriva- 
tion on how to expand the expectation value of a gen- 
eral operator as in Eq.([T]) in terms of density gradients. 
A well-known expansion method for a two-point func- 
tion is the Kirzhnits expansion 0] which is, however, 
specifically adapted to expanding the noninteracting one- 
particle density matrix in powers of density gradients and 
can not be used for more general correlation functions. 
The first goal of this paper is, therefore, to present a 
scheme based on nonlinear response theory that can be 
used to expand general correlation functions (such as the 
correlation hole) in terms of density gradients. 
The second goal of the paper is to clarify a point that 



is often overlooked in carrying out gradient expansions. 
We will illustrate the problem by considering the stan- 
dard gradient expansion for a global quantity, namely 
the exchange-correlation energy E xc . The starting point 
of any gradient expansion is the consideration of a den- 
sity profile n(r) = no + <5n(r) which consists of a small 
density variation 5n(r) around a homogeneous density 
uq. By considering the limit of slow density variations 
one then finds that the lowest gradient correction to the 
exchange-correlation energy is an integral over a func- 
tion of the form £? xc (no)(V<5n(r)) 2 where the coefficient 
B xc (no) is calculated from the static exchange-correlation 
kernel of the homogeneous electron gas of density no Q . 
Since V<5n(r) = Vn(r) we can replace the density vari- 
ation under the derivative operator by the full density 
profile n(r). However, the coefficient -B xc (no) still de- 
pends on the background density no. This is a problem 
in application of the formula to general inhomogeneous 
systems such as molecules or surfaces in which a back- 
ground density cannot be unambiguously defined (assum- 
ing that a low order gradient expansion makes sense in 
such systems [22|, HU). The usual argument to get rid 
of the background dependence, is that the replacement 
£? xc (no) — > B xc (no + <5n(r)) can be made since the differ- 
ence between these two quantities is or order Sn(r). How- 
ever, it is not clear that this is consistent with gradient 
coefficients that arise from terms of order (<5n(r)) 3 . The 
point was discussed clearly by Svendsen and von Barth 
[22| who checked that the replacement of no by the full 
density profile was consistent to third order in the density 
variation. This derivation was based on a specific rela- 
tion between response functions that describe the change 
in the exchange energy to second and third order in the 
density variations. In reference [24j this derivation was 
extended by showing that the replacement of no by the 
full density profile is consistent to all orders in the den- 
sity variation. In this paper we will show that this is 
the case as well for the gradient expansion of correlation 
functions rather than scalar functions. This relies on cer- 
tain relations between the response functions that must 
be satisfied for the gradient expansion to exist. 

The paper is divided as follows. In section |TT] we de- 
rive the basic equations and consistency conditions for 
the density gradient expansion of correlation functions. 
In section IIIII we derive the gradient expansion of the 
one-particle density matrix for a noninteracting system 
with density n(r) (which is therefore equal to the den- 
sity matrix of the Kohn-Sham system) and discuss its 
symmetry properties. We further calculate the gradient 
expansion of the exchange-hole to second order in the 
density gradients, both in real and in momentum space. 
The momentum space description is used to demonstrate 
that the system averaged exchange hole satisfies the sum 
rule (but not the negativity condition), and to calculate 
the gradient expansion of the exchange energy without 
the need to regularize divergent integrals. Finally in sec- 
tion [IV] we present our conclusions and outlook. 
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II. GRADIENT EXPANSION OF 
CORRELATION FUNCTIONS 

A. Expansion in density variations 

Let us for a system with ground state density n(r) con- 
sider an arbitrary correlation function F[n](r',r). Since 
the correlation function can be calculated from the many- 
body ground state, by virtue of the Hohenberg-Kohn the- 
orem this function is a functional of the density. At 
this point it will not be important what the specific 
form of the correlation function is, as the structure of 
their gradient expansions will be independent of its spe- 
cific form. Correlation functions of common interest are, 
for instance, the pair-correlation function, the exchange- 
correlation hole or the one-particle density matrix. 

When the density is constant in space, i.e. when 
n(r) = no, we are describing the homogeneous electron 
gas and due to translational and rotational invariance 
we have that -F[n](r',r) = F°(n ,|r — r'|) where F° 
is a function of the homogeneous density and the dis- 
tance between its spatial arguments. When the density 
n(r) = no + 5n(r) deviates from the constant density no 
by a small amount 5n(r) we can expand F in powers of 
this density variation. To derive compact expressions we 
first introduce the notation 

£m = ( r i i • ■ • > r m) 
dr m = rfri . . . dr m 

for the collection of m position vectors and the corre- 
sponding integration volume elements. We further define 

in 

3=1 

With these conventions the expansion of F in density 
variations is given by 

F[n](r',r) = F°(n , |r' - r|)+ 
00 1 f 

X) —\ \ dv m F m {n Q ;v' ,v,v m )5n{v m ) 

m=l ' ^ 

where we defined 



i^>o;r',r,r m ) 



(5 m F(r',r) 
<5n(ri) . ..Sn(r rn ) 



(4) 



The no dependence of the functions F m will be impor- 
tant for the gradient expansion, but to shorten notation 
we suppress this dependence in the notation for the time 
being and will reintroduce it later. Here we assume that 
the derivatives F m exist, which means that we assume 
that F is a smooth functional of the density. When the 



derivatives do exist we can expect the Taylor series (|4| 
to converge whenever |(5n(r)|/no <C 1. The density vari- 
ations are therefore assumed to be small but we do not 
need to put any constraint on how rapid they can vary in 
space and therefore their spatial derivatives can be large. 

Let us now look at the symmetry properties of the func- 
tions F m . As follows directly from the definition (j4} of 
the functions F m and the assumption of differentiability, 
we have the permutational symmetry 

F m (r', r, n . . . r m ) = F m (r\ r, r w(1) . . . r w(m) ) 

for all permutations ir of the integers 1 . . . m. Since the 
functions F m are evaluated at the homogeneous density 
no they also have the spatial symmetry properties of the 
homogeneous electron gas. These symmetries are the 
translational symmetry over an arbitrary vector a, rota- 
tional symmetry under arbitrary rotations by a rotation 
matrix i?, and inversion symmetry. If we denote by 

a = (a,..., a) 
Rr m = (Rri,...,Rr m ) 

the m-tuples of translation vectors and rotated position 
vectors we can write 

F m (r',r,rJ = F m (r> + a, r + a, r m + a) (5) 
^ m (r',r,r„J = F m (Rr',Rr,Rr m ) 
f™(r',r,r m ) = F m (-r', -r, -r m ). 

Since in Eq.([S]) the vector a is arbitrary we can in par- 
ticular choose a = r and define a new function N m 
depending on m + 1 independent vectors by the relation 



F m (r',r,n...r m ) 



F m (r'-r,0,r! -r, . 
iV m (r'-r,ri - r, . . 



r) 



(0) 



The difference vector r' — r will appear several times in 
our equations and it will therefore be convenient to define 
the short notation y = r' r. We will further use y = \y\ 
to denote the length of this vector. Our interest will 
be in the functions iV m and in particular their Fourier 
transforms 

^ m (y=qJ = J dv m N m {y,v_ m )e-^--^— (7) 



with Fourier inverse 



Ar m (n ;y,r m ) 



dq 



—m Arm 



(2tt) 3 



iqi-ri + ...+iq m -r„ 



With these definitions we obtain the following expansion 
for F in powers of the density variations 
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F[n](r',v) = F°(n , y) + >T ^ | dr m J — ^ N m (y t gJe^'^+--^<'^5n(T m ) 

00 1 f da 

= F°(n ,y)+ X; ~, J (^^ m (y.q m )e- i(qi+ - +q - ) - r ^(-q ro ) (8) 



This equation forms the basis of the gradient expansion. 
To proceed we further derive a consistency condition that 
is necessary for the existence of the gradient expansion. 
It follows directly from the definition (|4} of the functions 
F m that 

<5F m (r',r,r m ) = J dr" F m+1 (r', r, r", vJSn(v"). 

Taking 6n(r") = 6uq then to be a uniform shift in the 
density (which means that we look at a compression or 
decompression of the electron gas) yields 

f)F m f 

— (r',r,r m ) = J dr" F™+\r', r, r", rj. 

We can translate this condition to a condition on N m 
using its definition (|6|) and the definition (0 of its Fourier 
transform. This gives the condition 

8N m 

— — (y, qi . . . q m ) = A m+1 (y, 0, Ql . . . q ro ). (9) 
on 

This relation is of key importance for the existence of the 
gradient expansion. Without it we would not be able to 
eliminate the dependence of the gradient coefficients on 
the homogeneous density no in favor of a dependence on 
the actual density profile. As we will see, this requires 
a summation over response functions F m to infinite or- 
der. There is another important advantage of this re- 
summation. It will allow us to relax the constraint that 
the density variations be small (but varying arbitrarily 
fast) and replace it with the constraint that the density 
variations vary slowly but with an arbitrarily amplitude. 
This is discussed in the next section. 



B. The gradient expansion 

In order to expand F in density gradients we have to 
assume that the density gradients are small. This con- 
straint is most easily phrased in momentum space by re- 
quiring that the Fourier coefficients 6n(q) of the density 



variation have their main contribution from small wave 
vectors q. If this is the case then the main contribu- 
tion to the integrals in Eq.© comes from this region of 
small vectors and we can expand the functions N m in 
powers of the wave vectors. Subsequently carrying out 
the integrals in Eq. © then leads to an expansion in den- 
sity gradients, since powers in wave vectors correspond 
to orders of derivatives in real space. Since the response 
functions N m typically vary very rapidly for wave vec- 
tors close to the Fermi surface (or multiples thereof) we 
require that the Fourier coefficients Sn(q) of the density 
variation have their main contributions from wave vec- 
tors that satisfy |q| < kp where kp is the Fermi wave 
vector. 

The next task is then to expand the functions N m in 
powers of wave vectors. This task is simplified when we 
exploit the symmetry of the these functions. As follows 
directly from Eqs.© and ([7]) the functions TV™ in Fourier 
space inherit the symmetry properties of F m . They are 
symmetric functions of their wave vectors, i.e. 

N m (y, qi . . . q m ) = N m {y, . . . q, (m) ) 

for any permutation n of the integer labels, and are in- 
variant under rotations and inversion 

iV ro (y, qi ...q m ) = N m (Ry,Rq_ 1 ...R( lm ) 
7V™(y, qi ...q m ) = N m (-y,- qi ...-q m ). 

From these equations we therefore see that any power 
expansion of the functions N m must lead to an expansion 
in terms of the spatial vector y and the wave vectors 
qj which is invariant under rotations and inversions of 
these vectors. The mathematical question is then which 
polynomials have these properties. This was answered 
in the classic book by Weyl (2f|; every such polynomial 
must be a function of the variables y • y, y • q^ and q. ; ■ q^ . 
We see that these inner products are indeed invariant 
under rotations and inversions. Since any polynomial in 
powers of y 2 = y ■ y can be re-summed to a function of 
y we find that the general expression for N m to second 
order in the wave vectors q^ is given by 
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N m (y, q x . . q ro ) = N^(y) + JVTGO y • J> + iV 2 m (y) £q? + iV 3 m (y) £(y • q, 

?— 1 i— 1 6= 1 

rn rn 

+^ 4 m (y) • + ^r(y) £(y • *)(y ■*) + ••■ 



(10) 



i>j 



i>j 



where qf = ■ q^ and where wc took into account that 
this expansion must be invariant under permutations of 
the wave vectors q^. This expansion is readily extended 
to higher order powers in the wave vectors. For a given 
choice of correlation function F the practical task will be 
to determine the explicit form of the coefficient functions 
N" l (y) as a function of y and the background density no- 



If we insert Eq. (|10[) into Eq.® and Fourier transform 
back to real space we obtain the expansion 

OO -j oo 

F(r>, r) = F°(n , y) + £ — £ N?(y)A™(y, r). (11) 

m=l ' j=0 

Let us analyze the explicit form of the first five coeffi- 
cients A™ for j = 0, ..,4 of this expansion, since these 
are exactly the terms that we need for a gradient expan- 
sion to second order in the density derivatives. The first 
term for j = is simply given by 

A™ = Sn(r) m . 

The terms with j = 1, 2, 3 in Eg. pip involve according to 
Eq. (|10p the Fourier transforms of m symmetry equivalent 
terms and acquire the following forms in real space 



A? 
AT 



i m(Sn(r)) m 1 y • V<5n(r) 
-m(5n(r)) m - 1 V 2 c5n(r) 

-m(<5n(r)) m -V(- ■ V) 2 fo(r) 
V 



In the derivation of the last term we used that for an 
arbitrary function / 



■ Vj 



ViVj 



y y 



where with r = (xi,X2,xs) we used the notation di = 
d/dxi as well as yi = x[ — Xi. This expression is readily 
checked using 



ill 
y 



„3 



Finally the terms with j = 4, 5 in Eq.Qlip involve ac- 
cording to Eq. (|10p the Fourier transforms of m(m — l)/2 
symmetry equivalent terms which yields 



Af = --m(m 
1 



l)(<5n(r)) m - 2 (V5n(r)) 2 



AT = --m(m-l)(Sn(v)y 



- 2 (yVfe(r)f. 

We see that a general coefficient function A™ 1 depends 
on the density variation and gradients thereof. The 
functions Sn(r) that appear under the gradient opera- 
tors can be replaced by the actual density profile n(r) = 
no + 5n(r). However, we see that we are still left with 
powers of Sn(r) which are unprotected by derivative op- 
erators and which therefore can not be replaced by the 
full density profile n(r). It is precisely at this point that 
the consistency conditions (|9]) play a crucial role. If we 
use these conditions in Eq. (fT0| then we see that 



jym+1 

i 



(y) 



8NZ 



(»)■ 



(12) 



This equation relates certain coefficients coming from 
higher order response functions N m to those of lower 
order ones. In particular it tells us that by iteration 



N^(y) 



QmpO 

di 



iv) 



d m - l N] 

■^f(») 0-1.2,3) 
—^=f{y) (J=4,5). 



dn, 



If we insert these expressions together with the explicit 
form of the coefficient functions A™ back into Eq. (fTT]) in 
Eq- flT]) and shift some indices we obtain the expansion 
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°° 1 f) m 

F(r',r) = F\y) + £ ( y )5n(r)™ + 



m— 1 u 

> —5n(r) m \i— — 3-(y • Vn r ~ — -V n r) - — — 2_?/ 2 - • Vftir + 

m=0 o y 

-5 E [^< v "("' 2 + w (y - v " ,r))2) + (13) 

m=0 

In this expression we recognize several Taylor series of the coefficient functions A^ m (n +6n(r),y) in powers of 6n(r). 
These Taylor series can be re-summed to finally give 

F(r',r) = F°(n(r),y) + iNl(n(v),y)y Vn(r) - N^(n(r), y) V 2 n(r) - N^(n(r), y) y 2 (* ■ V) 2 n(r) 

y 

~N 2 (n(r), y) (Vn(r)) 2 - ^ 2 (n(r), y) (y • Vn(r)) 2 + . . . , (14) 

where the implicit dependence on no of the coefficient functions is now replaced by a dependence on the full density 
profile. We see that this is achieved by summing over response functions N m to infinite order. For clarity we reinstated 
the explicit density dependence of the coefficients in Eq. (fT4"l) to indicate that they are local functions of the density 
n(r) and therefore have a nontrivial spatial dependence on both r and y. We stress again that the derivative condition 
(|9]) was essential in eliminating the dependence of the gradient coefficient on the reference density no in favor of a 
dependence on the full density profile. 

Having obtained the general gradient expansion Eq. (|14p it remains to find explicit expressions for the coefficient 
functions NJ 1 . We will describe how to do this in detail for the coefficients needed for a gradient expansion to 
second order in the density derivatives. To calculate the coefficients N\, N 2 and N$ we need to calculate the function 
A'' 1 (no, y, q) for the homogeneous electron gas and expand it in powers of q: 

N\n , y, q) = <(n , y) + Nl(n , y)yq + N^(n , y)q 2 + N^(n , y)(y • q) 2 + . . . (15) 

The determination of the coefficients N 2 and iVf requires knowledge of the second order response function 

W 2 (n ,y, qi ,q 2 ) = N 2 (n , y) + N 2 (n , y)y ■ (qi + q 2 ) + N 2 (n , y)(q 2 + q 2 .) + 

jV 2 (n , y)((y • qi) 2 + (y • q 2 ) 2 ) + N 2 (n , y)( qi • q 2 ) + N 2 (n , y)(y • qi )(y • q 2 ) + . . . (16) 



The expression (|14j) together with Eqs. (fT5[) and (fT6|) de- 
termine completely the gradient expansion of an arbi- 
trary correlation function F to second order in the den- 
sity gradients. These equations form the main result of 
the present work. If expansions to higher order gradients 
are required they can be derived without difficulty along 
the lines described above. To do this one needs to con- 
struct higher order symmetric polynomials in the wave 
vectors and carry out the required Fourier transforms. 
The consistency conditions Eq.© or cquivalcntly Eq. (fT2)) 
then allow for a complete re-summation and leads to gra- 
dient coefficients depending on n(r). What is needed to 
determine the explicit form of the gradient coefficients in 
practice is the determination of the functions N m . How 
to do this for iV 1 and TV 2 is described in the next section. 



C. Determination of N 1 and N 2 

A practical calculation of the coefficients of the gradi- 
ent expansion of Eq. (fT4j) requires the knowledge of the 
the functions iV 1 and A^ 2 . According to Eqs.Q and © 
these are defined as the first and second density deriva- 
tives of the correlation function F with respect to the 
density, i.e. 



<SF(r',r) 
6n(r») lno 
SF(r',r) 
5n(r")8n{r'"y no 



iV 1 (y,r"-r) 
Ar 2 (y,r"-r,r"' - r). 



To derive useful expressions for these functions it is con- 
venient to transform derivatives with respect to the den- 
sity to derivatives with respect to the external potential 
as such derivatives appear naturally in perturbation the- 
ory. We can then use the well-developed tools of pertur- 
bation theory to calculate these functions. Let us there- 
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fore define the new functions 



THy r"-r) 

* ly ' r rj - Sv(r") lno 
S 2 F{r',r) 



_F 2 (y,r"-r,r'" - r) 



<5u(r")<M r '") 



(17) 



(18) 



Since we evaluate these functions for the homogeneous 
electron gas they only depend on differences between the 
spatial coordinates. It will therefore be convenient to also 
define their Fourier transforms by 



T (ri,r2) = J (2^^ (P ' q 
^^ r 3)=/^(k,p,q) 



ip-ri -Hq-ra 



ik-ri+ip-r 2 +iq-r 3 



as well as their partial Fourier transforms 

^ 1(y ' q)= /(l^ 1(P ' q)eiP ' y 



^ 2 (y,P,q)=/^ 2 (k,P,q) 



The chain rule of differentiation gives a connection be- 
tween the derivatives with respect to the density and the 
derivatives with respect to the potential of the correlation 
function F. We have 

SF(r u r 2 ) _ f ir ^f( ri ,r 2 ) Sn(v 4 ) 



5v(r 3 ) 
S 2 Ffa,r 2 ) 



5n(ri) Sv(r 3 ) 
5Ffa,r 2 ) <5 2 n(r 5 ) 



Svfa)Svfa) J ~~ Sn(r 5 ) 6vfa)5v(r 4 ) 
S 2 F(r 1 ,r 2 ) 6n(r 5 ) Snfa) 

dr 5 dr 6 



5n(r 5 )6n(r 6 ) Sv(r 3 ) <5u(r 4 ) ' 

We see that in these expressions the first and second order 
derivative of the density with respect to the potential 
appear. We therefore define the linear density response 
function by x 

Sn{ri '\ 
Si 



5n(ri) f dq 

mL = xfa - ri) = / — — jx(q 

)v(r 2 ) J (27r) d 



as well as the second order density response function \ 2 
by 

<5 2 n(ri) 2 

T-r-rn-^ ™ =x fa- ri,r 3 - n) 

dvfa)6vfa) 



(27T) 

Using these definitions we find by Fourier transforming 
Eqs.® and (gDJ) that 



i\T 1 (y q) = • Fl(y ' q) 
lY ' qj X(q) 



(21) 



^ 2 (y,P,q) 



[J 72 (y , P, q) - N 1 (y , p + q) X 2 (p, q)] 



x(p)x(q) 



(22) 



These equations give the desired relation between the 
density derivatives TV 1 and iV 2 and the potential deriva- 
tives T 1 and T 2 of the correlation function F . Relations 
between the higher order derivatives N m and F m can be 
derived in a completely analogous way. From Eqs. (|2"Tj) 
and (|22|) we see that to calculate the functions N 1 and 
N 2 , and hence the gradient expansion coefficients of F 
to second order in the gradients, we need to calculate 
the density response and the change in F to second 
order in a perturbing potential 8v. The problem is 
thereby reduced to doing perturbation theory on the 
electron gas. This is a well-developed field of research. 
One option is to use diagrammatic perturbation theory 
[26j |. In the following section we will use standard 
perturbation theory to obtain these response function 
for the case that F is the one-particle density matrix 
and use that to calculate the gradient expansion of the 
exchange hole. 



III. GRADIENT EXPANSION OF THE 
ONE-PARTICLE DENSITY MATRIX AND THE 
EXCHANGE HOLE 

A. The one-particle density matrix 

As an application of our formalism we carry out the 
gradient expansion of the one-particle density matrix for 
a noninteracting system with density n(r). This prob- 
lem has received large interest in the past since both the 
exchange energy E x [n] as well the Kohn-Sham kinetic 
energy T s [n] are explicit functionals of such a noninter- 
acting density matrix 0, Q • Therefore a gradient expan- 
sion of this density matrix directly leads to a gradient 
expansion of these two functionals. Such a gradient ex- 
pansion was first carried out by Gross and Dreizler [27j 
using the Kirzhnits expansion 0, HH . This expansion is 
adjusted to the specific form of the noninteracting one- 
particle density matrix and can not be generalized to ar- 
bitrary correlation functions. The method presented in 
this work can, on the other hand, be applied to arbitrary 
correlation functions, but to demonstrate the soundness 
of our formalism we will use it to provide an alternative 
derivation of the result obtained earlier by Gross and 
Dreizler using the Kirzhnits expansion. An advantage of 
our derivation based on nonlinear response theory is that 
it avoids the turning point problem encountered in the 
Kirzhnits approach [2J. 

Let us start by defining the one-particle density matrix 
in second quantization as 

7 (r>) = ^<^t( m )^( r V)|VI/) 

where a is a spin coordinate. We will consider the case 
of spin-compensated systems. Then, in a noninteracting 
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system with 27V electrons, the density matrix is given in 
terms of one-particle wave functions ifj (r) by 



N 

7(r',r) = ^2^(i>*(r) 

3=1 



(23) 



where the prc-factor 2 is a spin factor. The one-particle 
states are eigenstates of a one-particle Schrodingcr equa- 
tion 



~V 2 +v{r))<p j (r) = e j ip j (r) 



where v(r) is the external potential. To determine the 
gradient expansion coefficients of 7 we need to determine 
how the density matrix changes if we make small changes 
v(r) — > v(r) + 5v(r) in the external potential. More pre- 
cisely, we need to calculate the functional derivatives 



7 1 (ri,r 2 ,r 3 ) 



7 2 (ri,r 2 ,r 3 ,r 4 ) 



£7(1-1, r 2 ) 
Sv(r 3 ) 
S 2 j(r 1 ,r 2 ) 
8v(r 3 )6v(r4) 



which are the equivalents of the functions J- 1 and T 2 of 
Eqs. (fTT| and (TTS)) for the specific choice of correlation 
function F = 7 (note that they become functions of two 
and three independent vectors respectively when evalu- 
ated for a homogeneous system) . To do this it is sufficient 
to know the functional derivatives of the orbitals ipj and 
eigenvalues 6j with respect to the potential v. These are 
simply given by doing first order perturbation theory on 
the one-particle Schrodingcr Eq. (f24| . We find 



(25) 
(26) 



Sv(r) 

<^j( r ) 
Sv(r') 



l^(r)| 2 
' ^(r)^!-')^!-') 



E 



€j - e k 



With these relations we can differentiate Eq. (f23j) twice 
with respect to the potential. Afterwards we then need 
to insert the plane wave states an d their eigenval- 

ues ek appropriate for the homogeneous electron gas into 



the final expressions. In order not to interrupt the pre- 
sentation these calculations are given in Appendix IA1 and 
we simply present the results of these calculations here. 
If we define the Fermi factors by 

n k = 9(k F - |k|), 

where 8 is the Heaviside function and where &f = 
(37r 2 no) 1 / 3 is the Fermi wave vector, then the resulting 
expressions are given by 



7 (y>q) = 2 / 17. 

7 2 (y,p,q) = 2 J — e ik ^[$(k, k + p + q, k + p) 



(27) 



$(k,k 



where 



$(k,p,q) 



q,k 



(28) 



(29) 



In this expression e k = |k| 2 /2 are the one-particle ener- 
gies. 

It remains to calculate first and second order den- 
sity response functions x(q) an d X 2 (p, q) to evaluate the 
functions N 1 and N 2 of Eqs. flU) and fltfe]). Fortunately, 
for the specific choice F = 7 that we made these density 
response functions are simply given by 



x(q) = 7 1 (y = o J q) 
x 2 (p,q) = 7 1 (y = o,p,q). 



(30) 
(31) 



They are therefore obtained directly from Eqs. (|27p and 
(|2"5| so that we do not need to do any extra work to 
calculate them. 

To calculate the gradient coefficient to second order in 
the gradients we now need to expand 7 1 and 7 2 to second 
order in the wave vectors p and q. Since these calcula- 
tions are somewhat lengthy we do not present them here 
and refer to Appendix [B] for a description of the main 
steps. The resulting expansions are given by 



1 kp ( sin z 

7 (y,q) = --o 



7 2 (y,p,q) 



IT 2 \ Z 
1 



«q • y (q ■ y) : 



12k 2 



cos(z) 



2 

cos(z) - -(p + q) • y cos(z) + ^2 <- 



(p 2 + q 2 + p • q)(cos(z) + zsin(z)) 



24 



[(p • y) 2 + (q • y) 2 + 3((p + q) • y) 2 ] cos(z) 



(32) 



(33) 



where we defined z — kpy and q — |q|. The expansions for the first and second order response functions x an d X 2 
of Eqs. (f3"0")) and (|3T|) arc obtained by taking y = in these expressions. Together which Eqs. ([3"2")) and (|3"3")l these 
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expression then immediately determine the functions iV 1 and N 2 from Eqs. pTj) and (|22|) . The final result of this 
calculation to second order in the wave vectors is given by 



*q-y (q-y) 2 



^(y,P,q) = tt 



i'-'> 

F 



cos(z) — 



2 6 
sin(z) 



q 2 sin(z) — z cos(z) 
12k 2 z 



i-^(p + q)-y-^[(p-y) 2 + (q-y) 2 ] 



l 



12k 2 



(p 2 + q 2 + p • q)(3cos(z) + zsin(z) - 3 



sin(z) . 



Y^(p • y)(q • y)t 4 ^-^ - 3 «»(«)] 



By comparison of these two expression to the expansions ([15]) and ([16]) we can directly read off the gradient coefficient 
functions of Eq. (TT41 to be 



Nl(y) 



'do(z) Nl{y) 



12k 2 



12fc| 



[z 2 j {z)-Zz 3l (z)\ N 2 (y) 



12fc3 



[j {z) + 3zji(z)] 



where now we replaced tlq — > n(r) and hence the Fermi wave vector kp = (37r 2 n(r)) 1 / 3 that appears in these expressions 
is a spatially varying function. We further introduced the spherical Bessel functions 



JoO) 



sin z — z cos z 



By inserting the gradient coefficient function into Eq.(|f4j) we find the following explicit gradient expansion for the 
one-particle density matrix 



i 



i 



24fc5 



24k F 



(j (z)+3zj 1 (z))(y Vn(r)) 2 



(34) 



This expression is the main result of this section. After some manipulations we can rewrite it in an equivalent form 
as 



7(r',r) 



1 



1 1 



1 (Vfc 2 ) 2 
96 k% 



1 V 2 k 2 . . . 
24^T ZJ1(Z)+ I2fe 
f f 



32 k-p 



V ■ (Vfc| ■ ^) 

y 



(j (z)zJ - ZJl (z)) + -^r[Vk 2 ■ I ) (Z'jo(z) - Z^hiz)) 



y 2 ■ / n 

-z Jo{z) 



1. 



2 y 



zjo(z) 



(35) 



This expression becomes identical in form to that derived by Gross and Dreizler [27| using the Kirzhnits method 
after eliminating the effective Fermi vector of their work in favor of the density (this requires inversion of Eq.(f9) of 
reference [I?} )■ We close this section with a final remark on our result. We note that the gradient expansion to finite 
order breaks the symmetry 7(1*, r') = 7(1"', r)*. However, the full gradient expansion as described by Eq. ([TTj) has this 
symmetry which means that the symmetry is restored by taking all higher order gradients into account whenever the 
series converges. The symmetry violation has clearly consequences for the quantities that are derived from it, such as 
the exchange hole, as it introduces ambiguity in defining such functions. To be in accordance with existing literature 
we will in the next section adopt the definition of the exchange hole that was used by Perdew Q. 



B. The exchange hole and energy 

We finally stress a few points on the properties of the exchange hole as derived from the gradient expansion and 
present an alternative derivation of the gradient expansion for the exchange energy compared to that of Dreizler and 
Gross which has the advantage of avoiding the regularization of divergent integrals [27] . In doing so we confirm a 
result obtained by Langreth and Perdew The exchange hole can be calculated directly from the one-particle 

density matrix as 

| 7 (r + y,r)| 2 
^ (r + y|r) = My) ' 
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Inserting the expansion (|34[) of for the one-particle density matrix and retaining terms to second order in the gradients 
yields the expression 

Px(r + y|r) = -ln(r)iM _ 3 jo(*)ji(*) (y . Vn(r)) + _L Ji(z) 2 V 2 n(r) _ _L ZJo(z)jl(z)( ^ . V ) 2 n(r) 

+ ^M*K*Mz) ~ 3^(z))(Vn(r)) 2 + ^(^Mfl + 3. 7l (z) 2 - 3j (.) 2 )(y ■ Vn(r)) 2 + . . . (36) 

As was pointed out by Perdew [i| the second order gradient expansion of the exchange hole does not satisfy the 
negativity and sum rule constraints 

p x (r + y|r)<0 J dy p x (r + y|r) = -1. 

The violations of the negativity constraint is readily verified. However, the violation of the sum rule is not immediately 
obvious from Eq. (|36[) . The sum rule can, however, be conveniently analyzed in momentum space. To do this we define 
the Fourier transformed exchange hole to be 

p x (k|r)= J dy P Jr + y\r)e-*-y (37) 

This function has the form 

p x (k|r) =p°(k|r) + ai(fc, n) k • Vn(r) + a 2 {k, n) V 2 n(r) + a 3 (k, n) (k • Vn(r)) 2 

+ a 4 (fc, n)(V?i(r)) 2 + a 5 (k, n) (k • V) 2 ?i(r) + . . . (38) 

where we defined the unit vector k = k//c and k = |k|. The explicit form of the functions p x and aj follow directly b 
Fourier transforming Eq. (f3"6")l and are given in Appendix [Cj The coefficient functions aj are tempered distributions 
which have mathematically well-defined Fourier transforms. As follows directly from Eq. (|37[) the sum rule condition 
in momentum space translates to the requirement that /c x (k = 0|r) = —1. Since p x (k = 0|r) = —1 (sec Eq. (|Cl[) ) 
this implies that the sum rule would be fullfilled whenever <x, (fc = 0,n) = 0. However, as is clear from Eqs. (|C2[) - (|C6[) 
in the Appendix this is not satisfied for ai, whereas 013 and Q5 even diverge for k — > when interpreted locally as 
functions rather than distributions. The gradient expansion of the exchange hole does therefore indeed not satisfy the 
sum rule. These divergencies, however, cancel when we calculate the system averaged exchange hole in momentum 
space which is given by the expression 

N(p K (k)) = J drn(r) Px (k\r) 

= N(p°{k)) +k- [ drn(r) ai (k,n{r))S7n(r) + ^ f dr n(r) fa {k , n(r))a<n(r)a i n(r). (39) 
J i,j=i J 



where N is the number of electrons. In this expression functions aj as 
(p x (k)) is the system average of p x (k|r) and in the last 
term under the integral sign we defined the tensor 



CtT =0:4 



1 d(na 2 ) 



fejfc- / hk\ n dn 

fa(n, k) = a L (n, k)-^- + a T (n,k)[Sij % -f ) . 1 d(na 5 ) 

k z \ k z J a L =q t + "3 k 

n an 



This tensor has a longitudinal part with coefficient a^ 
and a transverse part with coefficient «t which describe 
the contributions to the system-averaged hole of density 
gradients Vn parallel and perpendicular to the momen- 
tum vector k. These coefficients are calculated from the 



as a short calculation on the basis of Eq. (f38|) will show. 
From these equations and the explicit expressions for aj 
given in Appendix [C] we can readily calculate the explicit 
expressions for q;l,t- If we define the dimcnsionless vari- 
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able k = k/(2kp) then we have 

"t =— , K Z T (k) 



QL 



32nfcf 



"32nfc| 



^L(fc), 



where the functions Zl,t have the explicit form 



Z T ( X ) = - 4x0(1 -x) + -S(x - 1) 



4x0(1 - x) + <J(s - 1) + ^«y'(a; - 1) 



These expressions are in accordance with the results of 
Langreth and Perdew (see Eq.(3.55) of reference [ljjj]) and 
this independent result therefore confirms the correctness 
of our alternative derivation. The interesting observation 
is now that 

lim Z LT (k) = 

fe^O 

and hence f3ij(k,n) — > when k — > 0. This implies that 
the last term in Eq.(f3"!)f does not contribute to the sum 
rule of the system averaged hole. The same conclusion 
can be derived for the second term in Eq. ([39|) . Since 
ai(k = 0,n) = z/(4nfcp) (see Eq. (|C2p ) is a local func- 
tion of the density the integrand of the first integral in 
Eg. ([59)) is a total derivative and vanishes (assuming ei- 
ther vanishing densities at infinity or periodic boundary 
conditions). We therefore find that 



lim(p x (k)) = lim(p x (k)) = 

fc->0 k— ¥0 



-1. 



This implies that system averaged exchange hole ob- 
tained from the second order gradient expansion does 
fulfill the sum rule, although the exchange hole itself docs 
not. We note, however, that this is only achieved by in- 
tegrating over both positive and negative contributions. 
When the positive contributions to the exchange hole are 
removed the sum rule is only recovered for a finite hole ra- 
dius [9|, [l0( . We finally note that with expression Eq. (|39|) 
it is now a simple matter to calculate the exchange energy 
from 



r , N 

EM = — 



dk 4tt 



We insert Eq. (j39|) and do the angular integrations first. 
It is therefore convenient to define the spherical and sys- 
tem averaged hole in momentum space by 



«ftc(*)» 



1 

47T 



dil k ( Px (k)) 



such that 



£*[»] = -/ dk((p x (k))). 



(40) 



From Eq. {[3j) we then find that 

N((p x (k))) = N((p x (k))) 



1 2 

drn(r)[-« L (fc,n) + -qt(^, n)](Vn) 2 



= N{{ P l(k))) 
where we defined 



dr 



Z x (k)(Vn) 2 



Z x (k) = ^(k) + ^Z T (k). 



(41) 



(42) 



The exchange energy is obtained by inserting the expres- 
sion for the averaged hole of Eq. (|4Tj) into Eq.(j40|). This 
then finally gives the expression 



E x [n] =- ^(^) 3 J rfrn3(r) 

+ [ dvB x {n{r)){Vn{r)) 2 



(43) 



where 



B x {n) 



16A;| Jo 



dk Z x (k) 
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4327r(37r 2 ) 1 /3 n 4 / 3 



and where to calculate the first term in Eq. (|4"3")l we used 
Eq. (|Cip . The coefficient B x is the same gradient coeffi- 
cient as obtained by Gross and Dreizler [27[ and earlier 
by Sham [3(J. However, the correct analytic exchange 
gradient coefficient is known (3lT[35| to be a factor 10/7 
larger. The reason for this discrepancy is clearly de- 
scribed by Svendsen and von Barth [34| who showed that 
the Coulomb interaction is too singular to allow for the 
interchange of the operations of integration and the ex- 
pansion in wave vectors. The problem does, for instance, 
not appear when Yukawa screened Coulomb interactions 
are used (3~H - The conclusion is that there is nothing 
wrong with the Kirzhnits method, or the nonlinear re- 
sponse theory derivation of the gradient expansion of the 
exchange hole presented here, but one has be aware that 
for Coulomb interactions the procedures of expanding 
correlation functions in terms of density gradients fol- 
lowed by integrations may not yield the same result as 
directly expanding the integrated quantity in terms of 
density gradients. For this reason the original GGA of 
Perdew and Wang [l(| based on the gradient ex pan sion 
of the exchange hole was later reparametrized jl6l . |36| 
to accomodate the correct gradient coefficient for the ex- 
change energy. 



IV. CONCLUSIONS AND OUTLOOK 

We derived a general scheme based on nonlinear re- 
sponse theory to calculate the density gradient expan- 
sion of general correlation functions and showed that in 
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order to express the gradient coefficients in terms of the 
full density profile summations to infinite order must be 
carried out over response functions of arbitrarily large or- 
der. A consistency condition was derived to do this. The 
formalism was used to derive the gradient expansion of 
the one-particle density matrix and the exchange hole to 
second order in the gradients. We confirm the derivation 
of Drcizler and Gross that used the Kirzhnits expansion 
method. We further analyzed the exchange hole in mo- 
mentum space to derive that the system averaged hole to 
second order in the gradients satisfies the sum rule and 
to derive the gradient expansion of the exchange energy 
without the need to regularize divergent integrals. 

The scheme that we presented is very general and can 
be applied to more general correlation functions. A nat- 
ural application of the scheme would be to calculate the 
gradient expansion of the correlation hole p c . Regarding 
the gradient terms of the correlation hole little is known. 
We essentially only have some detailed information on 
the long-range properties of the spherically and system 
averaged hole. This information comes from the work of 
Langreth and Pcrdew that calculated ((p c (k))) within the 
RPA. In our notation their results (see also Appendix C 
of [l9[) can be summarized as 

N(( Pc (k))) = N((p° c (k))) + J dr^z c (n,k)(Vnf 

where the function z c has been parametrized by Langreth 
and Mehl El EDl as 



z{n,k) 



(44) 



be desirable to have a first principles approach to the 
calculation of the short range properties of the gradient 
coefficient of the correlation hole. As follows from our 
derivation this requires the knowledge of the functions 
([TTJ) and (fT8f , or at least their expansion to second order 
in wave vectors, for the case that F represents the pair 
correlation function or the correlation hole of the electron 
gas. It may therefore be worthwhile to use our current 
scheme to explore these response functions beyond the 
RPA. For short range correlation an approach based on 
diagrammatic summation of ladder diagrams suggests it- 
self. Of course, for the development of density function- 
als for general systems beyond the weakly inhomogeneous 
regime it is not sufficient to use the straightforward gra- 
dient expansion. However, general short-range features, 
such as the way the exchange-correlation hole deforms 
close to the reference electron in the presence of a density 
gradient can be transferred to more general systems than 
the weakly inhomogeneous ones. Perhaps in the combi- 
nation with truly nonlocal functionals for the exchange- 
correlation hole [37} this can lead to the development of 
more accurate density functionals. This will be a topic 
of our future research. 
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where k s = {Ak-p/i^) 1 ' 2 is the Thomas-Fermi or screening 
wave vector. We can transform to real space to obtain 
the following expression for the spherically and system 
averaged correlation hole 



N{{p c (y)))=N{{p° c (y))) 
1 



dr 



964 (! 



(fc s2/ )Vl2) 2 



(Vn) 2 (45) 



We see from Eq. (j44)) that z c (n, k = 0) 7^ and as a con- 
sequence the sum rule ((p c (0))) = for the correlation 
hole is not satisfied. We know, however, that the RPA 
becomes accurate in the high density limit and since from 
Eq. (|45|) the gradient coefficient of the correlation hole de- 
pends on k s y we see that high densities are equivalent to 
large separations y. Similarly, low density corresponds 
to short-distance behavior. However, RPA can not be 
trusted in this region and we have no precise knowledge 
on the small distance behavior of the gradient coefficient 
of ((p c (y)))- A model for the general short-range behav- 
ior was proposed fltl - fl8| on physical arguments (it should 
not affect the Coulomb cusp of and the on-top value of 
the LDA hole) after which the real-space cutoff proce- 
dure was applied (see e.g. Fig. 4 in [17| and Fig. 5 in 
tl8f) to obtain a GGA for correlation. It would, however, 



Appendix A: Calculation of 7 1 and j 2 

In this Appendix we will derive in the expression for 
the first and second derivative of the one-particle density 
matrix 7 with respect to the potential v. For a clearer 
interpretation and compactness of notation we will only 
put a comma between the variables representing the orig- 
inal coordinates of the density matrix and the variables 
that appear as argument of the potential variations. Di- 
rect differentiation of Eq. (f2"B"| using Eq. (f2l))) then gives 



M23) 
Sv(l) 

OO 

2£(iW*) 



^■(1)^(1)^(2)^(3) 



(Al) 



where we used the short notation j = rj and introduced 
the occupation numbers fj = 1 for an occupied state and 
fj = for an unoccupied state. Inserting plane wave 
states appropriate for the electron gas we find that 

7 1 (r 2 r 3 ,ri) = 2 f rfqrfp np+q Z n P p ip-(r2-r 3 )-Hq-(ri-r 3 ) 



(27t)« 



tp+q 
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We therefore find that the function 7 1 has the simple 
form 



7 1 (p,q) = 2 



'p+q 



e p+q 



in Fourier space. This is precisely the integrand of 
Eq.pT]). The function 7 2 can be calculated by straight- 
forward differentiation of Eq. (|Aljl with respect to the po- 
tential 



7 2 (23,14) 



<V(23,1) 
dv(4) 



To obtain the explicit form of this function we have to dif- 
ferentiate both the orbitals and the eigenvalues using Eqs. 
(|25p and (f!?6| . Let us denote the term that arises from 
the differentiating the orbitals by X and the term that 
arises from differentiating the eigenvalues by Y . Then we 
find that 



7 2 (23, 14) = X(23, 14) + Y(23, 14), 



where 



(A2) 



X(23,U) = 2^^(2)^(3) 

$fe'fc)^(l)^(l)^(4)^(4) 
+ $Cj**)Vi(l)Vfc(lK(4)¥'*(4) 

F(23,14) = -2^^(2)^(3) 

$^)(|^-(4)| 2 -|^(4)| 2 )^(1)^(1) (A3) 

and where we further defined 



a- 



fk fi fj fi 



(A4) 



The function $ has the useful properties $>(ijj) = and 
®(ijk) = <&(ikj). The function 7 2 (23, 14) = 7 2 (23,41) is 
symmetric in the indices 4 and 1 due to the fact that the 
differentiations with respect to the potential commute. 
This is, however, not obvious from Eqs.(jA2| and (|A3[) . 
We therefore want to rewrite the form of the function 7 2 
in such a way that this symmetry is obvious. To do this 
we first note that since = the terms with 

k = i and k = j in Eq. (|A2[) sum up to a function Z of a 
similar form as Y in Eq. (|A3[) , namely 



2(23,14) = -2^> i (2) ¥ ,*(3) 



We can therefore write 7 2 (23, 14) as 

oo 

7 2 (23,14) = 2 J2 V<(2)v;(3) 

$fe-fcV fc (lK(lK(4)^-(4) 

f $0^)^(1)^(1)^(4)^(4) 
f F(23, 14) + Z(23, 14) 



(A6) 



It is easily seen that the sum of Y and Z is symmetric 
under interchange of 1 and 4. However, this is not obvi- 
ous in the first term of the equation above since at first 
sight $>(jik) does not appear to be equal to $(r/fc) since 



<&{jik) 



1 



fk - fj fi - fj 



Ck 



(A7) 



which seems to be different from expression (|A4p . How- 
ever, this is just appearance. The reason is that the oc- 
cupations can only attain the values zero and one. For 
the case fi = fj = 1 or /, = fj = we see directly that 
Eq. (|A4[) and (|A7|) attain the same value. A little inspec- 



tion shows that this is also true for cases /, = 0, fj = 1 
and /, = I, fj = 0. We therefore find for k ^ i,j that 
$>(ijk) = $(jifc). Therefore Eq. (|A6l) can be simplified to 



7 2 (23,14) = 2 ]T ^(2V*(3) 

ijk,k^(i,j) 

*(ijfc)[y*(l)^(l)^(4)^(4) 
+ ^■(1)^(1)^(4)^(4)" 



+ y(23, 14) + Z(23,14) 



(A8) 



This expression is now explicitly symmetric in the indices 
1 and 4. Let us now evaluate this expression for the 
homogeneous electron gas. In the electron gas the one- 
particle states are plane waves 

¥>k(r) = 



V 

where V is the volume of the system. Since |<^k| 2 = 1/V - 
it follows from Eqs. ([A3p and (|A5[) that the terms Y and Z 
are identically zero and the function 7(23, 14) of Eq. (|A8[) 
therefore attains the form 



7 2 (23,14) 



2 



$(k,p,q)e i ( k - r2 -P- r3 ) 

k,p,q,q#(k,p) 



e i(q-k)-ri+i(p-q)-r 4 , e 'i(p-q)-ri+i(q-k)-r 4 



where we defined 
$(k,p,q) = 

In this expression n F 



rik 



Ck 



fk 



'(ep — e p ) is the zero-temperature 

*(*J*)(l¥>j(l)l - IW 1 )! (4). (A5) Fermi function and £p = fc 2/ 2 ig the Fermi energy When 



14 



we replace the sum by an integral the restriction q^k,p 
is a region of measure zero and we can therefore freely 
integrate over all wave vectors. After a few substitutions 
we obtain 

7 * (23 , 14) = /^ 7 *( k , P , q > 

ik-(r 2 — r 3 )+ip-(ri— r 3 )+iq-(r4— r 3 ) 

where we defined 

7 2 (k,p,q) = 2[$(k,k+p+q,k+p)+$(k,k+p+q,k+q)] 
This gives the integrand of Eq. (|2"5f . 

Appendix B: Expansion of 7 1 and j 2 

1. Expansion of 7 1 

To determine 7 1 to second order in the wave vectors 
we have to expand the function 



(Bl) 



to second order in powers of q. This can be done by 
expanding the integrand in powers of q. If wc denote 
A = e p + q e p = p ■ q + q 2 /2 with q = |q|, then we 
can expand the Fermi function « p + q in a distributional 
Taylor series as 



dn, A 2 d 2 n 



A 3 dh 



n p+q = „ p ^ ^_ |ep ^ __ |ep ^ __ |ep ^ . . . 

Since n p = 9(e F — e P ) is the Heavisidc function we have 

A 2 . 

?i p+q =n p - A5(e F - e p ) + —5 (e F - e p ) 
A 3 

-—S"(e F -e p ) + ... 
Inserting this into Eq. (|Bl[) then gives 
7 1 (y,q) = -2/^3^ F -e p )e^ 

+ /(^ A5 ' (eF - £p)eiPy 

P r AV'(e F - ep )e^ + ... (B2) 



3 7 (2tt) 3 

We see that in these integrals several derivatives of the 
delta function appear. We now use the general mathe- 
matical expression 



where the sum runs over all zeros of the function y(x), i.e. 
y{xi) = 0. Using this equation in spherical coordinates 
with y(p) = (k F — p 2 )/2 we find 



S(e F 


-*p) 


S'(e F 


-e P ) 


6"(e F 


-e P ) 


S'"(e F 


-e P ) 



1 

k F 



1 



S'(p-k F ) 



pk F 

6 (p- k F ) 5" (p- k F ) 



k F p k F p 2 

1 r 1 "' 3 " 

—S (p-k F ) -5 (p - k F ) 



(B3) 
(B4) 

(B5) 



P 

\s'(p-k F ) 
p a 



P 



(B6) 



where in Eq. (|B6|) we also evaluated the third derivative of 
the delta function as we will need it later in the expansion 
of 7 2 . With Eqs. (|B3|) - (|B5j) we can readily evaluate the 
three integrals in Eq. (|B2[) . Let us denote these integrals 
by A, B and C. Then we find for the first integral 

k F J (2n)- i ir z z 

where z = k F y. The second integral gives 

B = j ^(p-q + 9 2 /2)<r(e F -e p )e^ 

= ( 9 2 /2 - *q • V y ) J j^S'(e F - e p ) 
1 



(q 2 /2 - iq ■ Vy) o _ 2) _ cos(fc F y) 
,,2 



47T 2 fcl 



cos z + 



yl 2n 2 k F 
iq-y sinz 



2tt 5 



-fcp- 



z 



where in the second step we used Eq. (|B4[) . It therefore 
remains to calculate the last integral 

C= -\j "(^(P^ + ^m.p-^e-- 
However, up to order q 2 it is suffices to calculate 

-(cos(z) + zsin(z)) 



67t 2 A;t 



■ cos(z) 



(q • y) 2 sin(z) 



67T 2 



k F - 



wherc in the second step wc used Eq. (|B5|) . Adding the 
results 7 1 = A + B + C gives the expression (f3"2"|) . 
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2. Expansion of 

Analogously to 7 1 we can expand 7 2 (y, p, q) of Eq. ([28| ) 
in powers of p and q. To reduce our effort we note that 
we can write 



7 2 (y, p, q) = ^(y, p, q) + ^(y, q, p) 



where 



f dk 

A(y,p,q)=2 J ^-^e lky $(k,k + p + q,k + p 



We therefore only need to expand A(y, p, q) in powers 
of the wave vectors and symmetrize with respect to p 
and q afterwards. To do this we first need to expand the 
function 

$(k,k + p + q,k + p) = 



1 



£k+p - £k+p+q 

If we denote 



"k+ P - Hk _ »k+ P +q - n k 

£k+p — £k Ck+p+q — £k 



Ai = e k + P - £k = — + P • k 
(p _)_ q)2 

A 2 = ek+p+q - e k = ^ h ( p + q ) ■ k ' 

we find by expanding the Fermi functions in a distribu- 
tional Taylor series that 

$(k, k + p + q, k + p) = -S' (e F - e k ) 

- -(Ai + A 2 )S"(e F -e k ) 
6 

+ ^(A 2 + AiA 2 + A^"'(e F - e k ) + . . . 



With this expression we find that the function A(y, p, q) 
can be calculated as 



A(y,p,q) 

4/(^ y ' (eF ^ k)(Ai+A2)e " y 

+ hj (^3 -^(Aj + AiAa + ADe^ 



The three integrals appearing on the left hand side of this 
expression are sufficient to extract the powers to second 
order in p and q. Let us call these three integrals Ax, A% 
and A3 in order of appearance. The first integral gives 



Ax 



dk 

(2^ 



S (e F - e k ) e lk - 



iky 



27T 2 fc F ' 



The second integral is given by 



dk 



(2tt) 3 

P 2 + (P + q) 2 
6 

p 2 + (p + q) 2 

6 



S (e F - £k) 



p 2 + (p + q) 5 



(2p + q) • k 



-(2p + q)-V y 
^(2p + q)-V y 



cos(z) + z sin(z) 



-2tt 2 /c f 



cos(z) + zsin(z) b2 + (p + q)2) 



127T 2 fc| 



(2p + q)-y cos(z) 



67T 2 fc F 

where we used Eq. (|B5P in the second step. Finally the third integral has the explicit form 



-4, 



1 



dk 



+ ( 



12 J (2tt) 3 

(p + q) 2 , 



4 + P-k) 2 + 4 + p.k)(iE±S)! 

(p + q)-k) 2 lr(e F - ek )e lk -y 



+ (p + q) • k) 
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However, we only need the terms to second order in p and q. Therefore, it is sufficient to calculate 
1 P dk 

Ai = 12 J j^r [(p ' k)2 + (p ' k)((p + q) ■ k) + ((p + q) ■ k)2] 6 {ep £k) elky 

1 r dk 

= [(P ' v y) 2 + (P ' v y)((P + q) ' v y) + ((p + <0 ' v y) 2 ] / J^s 6 ^ " e ^ e y 

1 r . _ , 9 , _ ... , ^ j . ,. . _ . 91 /3cos(z) — z 2 cos(z) + 3zsin(z)\ 

= [(P ' V y ) 2 + (p ■ V y )((p + q) ■ V y ) + ((p + q) ■ V y ) 2 ] ^ ^ 



where to evaluate the integral over the delta function we used Eq. (|B6p . To perform the derivatives in the last term 
we use that for an arbitrary function f(z) 

~d 2 f Idf 



(p.V y )(q.V y )/(.) = ( P .y)(q.y)4 

y 



dz 2 z dz 



+ (p-<l) — -r- 

V dz 



With this expression we find that 

1 

24?T 2 fc F 
1 



A 3 = - WJZvU (P 2 + P ■ (P + q) + (P + q) 2 )(cos(^) + zsm(z)) 

[(p • y) 2 + (p • y)((p + q) ■ y) + ((p + q) ■ y) 2 ] cos(z) 



247T 2 fc F ' 

Collecting our results A = A\ + A2 + A3 we therefore obtain the following expression for 7 2 (y, p, q) 

7 2 (y,P,q) = A(y,p,q)+A(y,q,p) 

1 / i 1 

= ( cos(z) - -(p + q) • ycos(z) + ^^"(P 2 + q 2 + P ' q)(cos(z) + zsin(z)) 

-^[(p • y) 2 + (q • y) 2 + 3((p + q) • y) 2 ] cos(z) 



Appendix C: Exchange hole in momentum space 



This is equal to expression (|33 



Here we present the explicit expressions for the coefficient functions aij of Eq. (|38| which are calculated by Fourier 
transforming Eq. (|36|) . If we define the dimcnsionless quantity k = |k|/(2fcp(r)) then p° is given by 

^(k|r) = (-l + ^fc-ifc 3 Ml-fc) (CI) 

which is simply the Fourier transform of the LDA hole. We see that p x (0|r) = —1 and therefore the LDA hole satisfies 
the sum rule. The functions ctj are given by 

ax=^(l-*) (C2) 



^=-^0il-k) (03) 



F 



*--H^( ? T^ + i'< t - 1 > + i'< t - 1 >) (C4) 



"F 
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We present a general scheme based on nonlinear response theory to calculate the expansion of 
correlation functions such as the pair-correlation function or the exchange-correlation hole of an 
inhomogeneous many-particle system in terms of density derivatives of arbitrary order. We further 
derive a consistency condition that is necessary for the existence of the gradient expansion. This 
condition is used to carry out an infinite summation of terms involving response functions up to 
infinite order from which it follows that the coefficient functions of the gradient expansion can be 
expressed in terms the local density profile rather than the background density around which the 
expansion is carried out. We apply the method to the calculation of the gradient expansion of the 
one-particle density matrix to second order in the density gradients and recover in an alternative 
manner the result of Gross and Dreizler (Z. Phys. A 302, 103 (1981)) which was derived using 
the Kirzhnits method. The nonlinear response method is more general and avoids the turning 
point problem of the Kirzhnits expansion. We further give a description of the exchange hole in 
momentum space and confirm the wave vector analysis of Langreth and Perdew (Phys. Rev. B21, 
5469 (1980)) for this case. This is used to derive that the second order gradient expansion of the 
system averaged exchange hole satisfies the hole sum rule and to calculate the gradient coefficient 
of the exchange energy without the need to regularize divergent integrals. 
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I. INTRODUCTION 



Since the groundbreaking work of Hohenberg and 
Kohn [l[ we know that the external potential of any inho- 
mogeneous quantum many-body system is a functional of 
its ground state density n. This implies that the many- 
body ground state and hence any ground state 
expectation value 



as 



0,a 



A[n] = (*[n]|i|*[n]) 



(1) 



for any operator A is a functional of the density. This 
idea has inspired an enormous amount of work in a re- 
search field that is now known as density-functional the- 
ory. Density-functional theory 041] is currently one of 
the main approaches used in electronic structure the- 
ory. Over the years many extensions beyond its origi- 
nal formulation have been developed and currently it is 
widely applied in solid state physics and quantum chem- 
istry, both for ground state and time-dependent systems 
Q. Especially a large activity has gone into finding ex- 
plicit expressions for the ground state energy E[n] as a 
functional of the density for many-electron systems. The 
ground state energy is obtained by minimizing this func- 
tional on an appropriate set of electronic densities, which 
is of great importance in determining structures and ge- 
ometries of molecules and solids. By use of the Kohn- 
Sham method @ this minimization problem is reduced 
to solving effective one-particle equations. The construc- 
tion of the effective potential in these equations requires 
the knowledge of the so-called exchange-correlation en- 
ergy functional E xc [n] . This functional can be expressed 



E*M = 5 J dvdv | r _ r ,| 



(2) 



where p xc (r'|r) is the coupling constant integrated 
exchange-correlation hole. The exchange-correlation hole 
has a physical interpretation as the difference between 
the conditional and the unconditional probability (which 
is simply the density n(r')) to find a electron at r', given 
that we know that there is an electron at r. Equation ([2|) 
is an important relation since it expresses the exchange- 
correlation energy in terms of a quantity that has a lo- 
cal physical interpretation and can be studied by accu- 
rate wave function and many-body methods or modeled 
based on physical intuition. It has therefore played an 
important role in the development of approximate den- 
sity functionals (cj-[T^|. One of the simplest and already 
very successful approximations has been the local density 
approximation (LDA) in which the exchange-correlation 
hole is taken from the homogeneous electron gas and ap- 
plied locally [|| to systems with arbitrary density profiles. 
The accuracy of the LDA has been considerably improved 
by means of the so-called generalized gradient approxi- 
mations (GGAs) [l5|. A large class of commonly used 
GGAs is based on the so-called real-space cutoff of the 
straightforward gradient expansion for the exchange and 
correlation hole. The first such functional for exchange 
was proposed by Perdew Q . He noted that the gradient 
expansion of the exchange hole to second order in the 
density gradients violates the negativity condition and 
the hole sum rule by which the exchange hole integrates 
to a deficit of one electron. By enforcing these physical 
constraints a density functional was obtained that greatly 



2 



improved on the LDA for binding energies of molecules. 
This procedure was later simplified [lfj using the fact 
that the exchange correlation energy ([2]) can be written 

as 



N f 
N = y J dy 



(pxc(y)) 



where N is the number of electrons in the system and 
(p X c(y)) = Jj J dr?i(r)p xc (r + y|r), (3) 

is the so-called system-averaged exchange-correlation 
hole. This averaged hole still obeys the negativity con- 
dition as well as the sum rule and can therefore be sub- 
jected to the real-space cut-off procedure. One advan- 
tage of using the system averaged holes was to reduce 
the order of the derivatives in the gradient expansion. 
Furthermore it also allowed for the real-space cut-off pro- 
cedure to be applied to the correlation hole since there is 
a known gradient expansion for the system averaged cor- 
relation hole in the random phase approximation (RPA) 
but no known expansion for the correlation hole itself. 
This was the basis for a GGA for the correlation energy 
SEMI- These GGAs relied heavily on the pioneering 
work by Langreth, Perdew and Mehl [l9l - |2~lj ] who per- 
formed a wave-vector decomposition of the system aver- 
aged hole of Eq. and calculated the Fourier transform 



(Pxc(k)) = / dy(p xc {y))e 



-iky 



from which the exchange correlation energy can be cal- 
culated as 



N 



We are not aware of any work that goes beyond refer- 
ences [l9l - |2~l| and improves on the straightforward gra- 
dient expansion of the system averaged correlation hole. 
The knowledge on the gradient expansion of the correla- 
tion hole is very limited indeed. As mentioned before, in 
contrast to the work on the exchange hole 0, @, there 
are no known expressions for the gradient expansion of 
the non-system averaged correlation hole. 
Part of the problem is that there has been no clear deriva- 
tion on how to expand the expectation value of a gen- 
eral operator as in Eq.([T]) in terms of density gradients. 
A well-known expansion method for a two-point func- 
tion is the Kirzhnits expansion 0] which is, however, 
specifically adapted to expanding the noninteracting one- 
particle density matrix in powers of density gradients and 
can not be used for more general correlation functions. 
The first goal of this paper is, therefore, to present a 
scheme based on nonlinear response theory that can be 
used to expand general correlation functions (such as the 
correlation hole) in terms of density gradients. 
The second goal of the paper is to clarify a point that 



is often overlooked in carrying out gradient expansions. 
We will illustrate the problem by considering the stan- 
dard gradient expansion for a global quantity, namely 
the exchange-correlation energy E xc . The starting point 
of any gradient expansion is the consideration of a den- 
sity profile n(r) = no + <5n(r) which consists of a small 
density variation 5n(r) around a homogeneous density 
uq. By considering the limit of slow density variations 
one then finds that the lowest gradient correction to the 
exchange-correlation energy is an integral over a func- 
tion of the form £? xc (no)(V<5n(r)) 2 where the coefficient 
B xc (no) is calculated from the static exchange-correlation 
kernel of the homogeneous electron gas of density no Q . 
Since V<5n(r) = Vn(r) we can replace the density vari- 
ation under the derivative operator by the full density 
profile n(r). However, the coefficient -B xc (no) still de- 
pends on the background density no. This is a problem 
in application of the formula to general inhomogeneous 
systems such as molecules or surfaces in which a back- 
ground density cannot be unambiguously defined (assum- 
ing that a low order gradient expansion makes sense in 
such systems [22|, HU). The usual argument to get rid 
of the background dependence, is that the replacement 
£? xc (no) — > B xc (no + <5n(r)) can be made since the differ- 
ence between these two quantities is or order Sn(r). How- 
ever, it is not clear that this is consistent with gradient 
coefficients that arise from terms of order (<5n(r)) 3 . The 
point was discussed clearly by Svendsen and von Barth 
[22| who checked that the replacement of no by the full 
density profile was consistent to third order in the density 
variation. This derivation was based on a specific rela- 
tion between response functions that describe the change 
in the exchange energy to second and third order in the 
density variations. In reference [24j this derivation was 
extended by showing that the replacement of no by the 
full density profile is consistent to all orders in the den- 
sity variation. In this paper we will show that this is 
the case as well for the gradient expansion of correlation 
functions rather than scalar functions. This relies on cer- 
tain relations between the response functions that must 
be satisfied for the gradient expansion to exist. 

The paper is divided as follows. In section |TT] we de- 
rive the basic equations and consistency conditions for 
the density gradient expansion of correlation functions. 
In section IIIII we derive the gradient expansion of the 
one-particle density matrix for a noninteracting system 
with density n(r) (which is therefore equal to the den- 
sity matrix of the Kohn-Sham system) and discuss its 
symmetry properties. We further calculate the gradient 
expansion of the exchange-hole to second order in the 
density gradients, both in real and in momentum space. 
The momentum space description is used to demonstrate 
that the system averaged exchange hole satisfies the sum 
rule (but not the negativity condition), and to calculate 
the gradient expansion of the exchange energy without 
the need to regularize divergent integrals. Finally in sec- 
tion [IV] we present our conclusions and outlook. 
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II. GRADIENT EXPANSION OF 
CORRELATION FUNCTIONS 

A. Expansion in density variations 

Let us for a system with ground state density n(r) con- 
sider an arbitrary correlation function F[n](r' ,r). Since 
the correlation function can be calculated from the many- 
body ground state, by virtue of the Hohenberg-Kohn the- 
orem this function is a functional of the density. At this 
point it will not be important what the specific form of 
the correlation function is, as the structure of its gradient 
expansion will be independent of its specific form. Corre- 
lation functions of common interest are, for instance, the 
pair-correlation function, the exchange-correlation hole 
or the one-particle density matrix. 

When the density is constant in space, i.e. when 
n(r) = no, we are describing the homogeneous electron 
gas and due to translational and rotational invariance 
we have that -F[n](r',r) = F°(n ,|r — r'|) where F° 
is a function of the homogeneous density and the dis- 
tance between its spatial arguments. When the density 
n(r) = no + 5n(r) deviates from the constant density no 
by a small amount 5n(r) we can expand F in powers of 
this density variation. To derive compact expressions we 
first introduce the notation 



r 

— n 



(ri,...,r r 
dri . . . dr Tl 







for the collection of m position vectors and the corre- 
sponding integration volume elements. We further define 



With these conventions the expansion of F in density 
variations is given by 

F[n](r',r)=F (n 0j |r'-r|)+ 



E 

m— 1 



dr m F m (n Q ;r',r,r m )5n(r m ) 



where we defined 



F™(n ;r',r,r m ) 



S m F(r', r) 
5n(ri) . ..Sn(r m ) 



(4) 



The n dependence of the functions F m will be important 
for the gradient expansion, but to shorten notation we 
suppress this dependence in the notation for the time be- 
ing and will reintroduce it later. Here we assume that the 
derivatives F m exist, which means that we assume that 
F is a smooth functional of the density. In the electron 
gas two Taylor expansions around two different values of 
no have very likely the same mathematical convergence 
properties since the two systems have the same physical 
properties (unless we are close to a very low density cor- 
responding to a Wigner crystal transition) . The point uq 



is therefore unlikely to be a non- analytic point. When 
the derivatives F m do exist we can therefore expect the 
Taylor series (j4} to converge whenever \5n(r)\/n < 1 or 
when they are integrable with a small integral norm. The 
density variations are therefore assumed to be small but 
we do not need to put any constraint on how rapid they 
can vary in space and therefore their spatial derivatives 
can be large. We note, however, that the notion of the 
existence of the functional derivative is closely related 
to the allowed domain of densities and the norm defined 
on this function space. The smallness or integrability of 
|<5n(r)| implies the use of a supremum norm or possibly 
an i p -integral norm (i.e. requiring |<5n(r)| p to be inte- 
grable). Other norms may induce weaker constraints on 
the derivative functions F m but stronger constraints on 
the density variations Sn(r). A rigorous discussion on 
the issue of functional differentiability for the case of the 
energy functional is given in Refs. [25l - l28j . 

Let us now look at the symmetry properties of the func- 
tions F m . As follows directly from the definition Q of 
the functions F m and the assumption of differentiability, 
we have the permutational symmetry 

F m (r', r, n . . . r m ) = F m (r\ r, r w(1) . . . r w(m) ) 

for all permutations ir of the integers 1 . . . m. Since the 
functions F m are evaluated at the homogeneous density 
no they also have the spatial symmetry properties of the 
homogeneous electron gas. These symmetries are the 
translational symmetry over an arbitrary vector a, rota- 
tional symmetry under arbitrary rotations by a rotation 
matrix R, and inversion symmetry. If we denote by 

a = (a,..., a) 
Rr m = (Rri,...,Rr m ) 

the m-tuples of translation vectors and rotated position 
vectors we can write 



F m (r',r,rJ = F m (r' + a, r + a, r r 
F m (v',r,v m ) = F m (Rr\ Rr, RrJ 
f™(r',r,r m ) = F m (-r', -r, -rj. 



(5) 



Since in Eq.([S]) the vector a is arbitrary we can in par- 
ticular choose a = r and define a new function N m 
depending on m + 1 independent vectors by the relation 



F m (r',r,n...r m ) 



F"V-r,0,n-r,...,r m -r) 
JV m (r'-r,n -r,...,r m -r). 



(6) 



The difference vector r' — r will appear several times in 
our equations and it will therefore be convenient to define 
the short notation y = r' r. We will further use y = |y| 
to denote the length of this vector. Our interest will 
be in the functions N m and in particular their Fourier 
transforms 



*Hy.aJ= / * m iV™(y,r m )e- 



jqi-ri- 



-iq„ 



(7) 
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with Fourier inverse With these definitions we obtain the following expansion 

for F in powers of the density variations 

iV-(n ;y,r m ) = J N ™(y, qje— . 



F[n](r',r) = FVo,2/)+£^y / *m / q m )e^-( r - r )+-+^-^- r )«5n(r m ) 



m=l 
oo 



00 1 /• <iq 



(8) 



This equation forms the basis of the gradient expansion. 
To proceed we further derive a consistency condition that 
is necessary for the existence of the gradient expansion. 
It follows directly from the definition (|4]) of the functions 
F m that 

SF m (v',v,v m ) = J dr" F m+ \r> ,r,r" ,r m )Sn(v"). 

Taking 5n(r") = 5no then to be a uniform shift in the 
density (which means that we look at a compression or 
decompression of the electron gas) yields 

— {r',r,r_ m ) = J dr" F™+\r' , r, r" , rj. 

We can translate this condition to a condition on N rn 
using its definition ([6]) and the definition (0 of its Fourier 
transform. This gives the condition 

8N m 

— (y, q x . . . q m ) - A m+1 (y, 0, qi . . . q m ). (9) 

dn 

This relation is of key importance for the existence of the 
gradient expansion. Without it we would not be able to 
eliminate the dependence of the gradient coefficients on 
the homogeneous density n in favor of a dependence on 
the actual density profile. As we will see, this requires 
a summation over response functions F m to infinite or- 
der. There is another important advantage of this re- 
summation. It will allow us to relax the constraint that 
the density variations be small (but varying arbitrarily 
fast) and replace it with the constraint that the density 
variations vary slowly but with an arbitrary amplitude. 
This is discussed in the next section. 



B. The gradient expansion 

In order to expand F in density gradients we have to 
assume that the density gradients are small. This con- 
straint is most easily phrased in momentum space by re- 
quiring that the Fourier coefficients Sn(q) of the density 
variation have their main contribution from small wave 



vectors q. If this is the case then the main contribu- 
tion to the integrals in Eq.© comes from this region of 
small vectors and we can expand the functions N m in 
powers of the wave vectors. Subsequently carrying out 
the integrals in Eq. © then leads to an expansion in den- 
sity gradients, since powers in wave vectors correspond 
to orders of derivatives in real space. Since the response 
functions N m typically vary very rapidly for wave vec- 
tors close to the Fermi surface (or multiples thereof) we 
require that the Fourier coefficients Sn(q) of the density 
variation have their main contributions from wave vec- 
tors that satisfy |q| < k-p where k-p is the Fermi wave 
vector. Note that the procedure of interchanging inte- 
gration and summation requires absolute convergence of 
the wave vector expansion. It is hard to prove this prop- 
erty for the general response functions that we consider 
here and we therefore have to assume its validity. 

The next task is then to expand the functions N m in 
powers of wave vectors. This task is simplified when we 
exploit the symmetry of the these functions. As follows 
directly from Eqs.© and ([7]) the functions N m in Fourier 
space inherit the symmetry properties of F m . They are 
symmetric functions of their wave vectors, i.e. 

N m (y, qi . . . q m ) = N m (y, q ff(1) ...q ff ( m )) 

for any permutation n of the integer labels, and are in- 
variant under rotations and inversion 

A m (y, qi ...q m ) = N m (Ry,R qi ...Rq m ) 
A m (y, qi ...q m ) = AH-y,-qi...-q ro ). 

From these equations we therefore see that any power 
expansion of the functions A" 1 must lead to an expansion 
in terms of the spatial vector y and the wave vectors 
qj which is invariant under rotations and inversions of 
these vectors. The mathematical question is then which 
polynomials have these properties. This was answered 
in the classic book by Weyl (29|; every such polynomial 
must be a function of the variables y • y, y • q^ and q^ • q^ . 
We see that these inner products are indeed invariant 
under rotations and inversions. Since any polynomial in 
powers of y 2 = y • y can be re-summed to a function of 
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y we find that the general expression for N m to second order in the wave vectors is given by 



J 



7V m (y, q x . . . q m ) = N^(y) + N{ n (y) y • »< £ + X> ' * 

6=1 i—1 i—1 

rn rn 

+N?(y) £(<b • Cfc) + iV™(y) £(y • q ? )(y ■*) + ... 



(10) 



t>3 



r 



where qf = qi ■ qi and where we took into account that 
this expansion must be invariant under permutations of 
the wave vectors q^. This expansion is readily extended 
to higher order powers in the wave vectors. For a given 
choice of correlation function F the practical task will be 
to determine the explicit form of the coefficient functions 
N" l (y) as a function of y and the background density no- 



If we insert Eq. (|10[) into Eq.([5]) and Fourier transform 
back to real space we obtain the expansion 



OO -. oo 

F(r>, r) = F°(n , ?/)+£-£ N?(y)A™(y, r). (11) 



V 

^ mi 

m=l j=0 



Let us analyze the explicit form of the first six coefficients 
A™ 1 for j = 0, .., 5 of this expansion, since these are ex- 
actly the terms that we need for a gradient expansion to 
second order in the density derivatives. The first term 
for j = is simply given by 



Sn(r) r 



The terms with j = 1, 2, 3 in Eq. ([TT]) involve according to 
Eq. (fTU| the Fourier transforms of m symmetry equivalent 
terms and acquire the following forms in real space 



A™ 
A™ 



i m(Sn(r)) m 1 y • V<5n(r) 
-m((5n(r)) m - 1 V 2 (5n(r) 

-m(5n(r)) m "V(- ■ V) 2 5n(r) 
V 



In the derivation of the last term we used that for an 
arbitrary function / 



■ Vj 



ViVj 



y y 



v- 



where with r = (x\,X2,xs) we used the notation di = 
d/dxi as well as yi = x[ — xi. This expression is readily 
checked using 



ill 
y 



„3 



Finally the terms with j = 4, 5 in Eq. (|ll[) involve ac- 
cording to Eq. (|10[) the Fourier transforms of m(m — l)/2 
symmetry equivalent terms which yields 



-4', 



A' 



m(m - l)(<5n(r))"^ 2 (VSn(r)) 



m(m - l)(<5n(r)) m ^ (y • V5n(r)) 



We see that a general coefficient function AJ 1 depends 
on the density variation and gradients thereof. The 
functions Sn(r) that appear under the gradient opera- 
tors can be replaced by the actual density profile n(r) = 
no + Sn(r). However, we see that we are still left with 
powers of Sn(r) which are unprotected by derivative op- 
erators and which therefore can not be replaced by the 
full density profile n(r). It is precisely at this point that 
the consistency conditions (|9]) play a crucial role. If we 
use these conditions in Eq. ([T0| then we see that 



' J 



(y) 



dNV 
dn 



(y)- 



(12) 



This equation relates certain coefficients coming from 
higher order response functions N m to those of lower 
order ones. In particular it tells us that by iteration 



N^iy) 



QmpO 



dnl 



(y) 



d m - l N} 

-^-(y) = 1,2,3) 

-rs=f(l/) (7 = 4,5). 



dn 



If we insert these expressions together with the explicit 
form of the coefficient functions AJ 1 back into Eq. (fTTj) in 
Eq. flT]) and shift some indices we obtain the expansion 
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°° 1 f) m 

F(r',r) = F\y) + £ ( y )5n(r)™ + 



m—l u 



oo 



S ^ ?l(r) ^W (y ' W(r)) " V n(r) " { v ' V) n(r)] + 

-5 E s 4 "'"™ [^< v "("' 2 + w (y - v " ,r))2) + (13) 

m=0 

In this expression we recognize several Taylor series of the coefficient functions Np(n +6n(r),y) in powers of 6n(r). 
These Taylor series can be re-summed to finally give 

F(r',r) = F°(n(r),y) + iNl(n(v),y)y Vn(r) - N^(n(r), y) V 2 n(r) - (n(r) , y) y 2 ■ V) 2 n(r) 

y 

~N 2 (n(r), y) (Vn(r)) 2 - ^ 2 (n(r), y) (y • Vn(r)) 2 + . . . , (14) 

where the implicit dependence on uq of the coefficient functions is now replaced by a dependence on the full density 
profile. We see that this is achieved by summing over response functions N m to infinite order. For clarity we reinstated 
the explicit density dependence of the coefficients in Eq. (fT4")) to indicate that they are local functions of the density 
n(r) and therefore have a nontrivial spatial dependence on both r and y. We stress again that the derivative condition 
(|9]) was essential in eliminating the dependence of the gradient coefficient on the reference density no in favor of a 
dependence on the full density profile. Having obtained the general form of gradient expansion we can wonder about 
its convergence. In order for the gradient expansion to be useful it at least needs to be an asymptotic series. This 
is known to be the case for the gradient expansion of some commonly studied functionals. For instance, for small 
densities variations very accurate results for the exchange energy are obtained by summing all terms up to fourth 
order in the density derivatives [3, HH . 

Having obtained the general gradient expansion Eq. (|14l) it remains to find explicit expressions for the coefficient 
functions N™. We will describe how to do this in detail for the coefficients needed for a gradient expansion to 
second order in the density derivatives. To calculate the coefficients N^N^ and N% we need to calculate the function 
iV^no, y, q) for the homogeneous electron gas and expand it in powers of q: 

JV^no, y, q) = iVoOo, y) + Nl{n , y)yq + N^(n , y)q 2 + JV^no, y)(y • q) 2 + • • • (15) 

The determination of the coefficients N 2 and N 2 requires knowledge of the second order response function 

7V 2 (n ,y, qi ,q 2 ) = N 2 (n , y) + N 2 (n , y)y • (q a + q 2 ) + N 2 (n a , y)(q 2 + q 2 ) + 

iV 2 (n , y)((y • qi) 2 + (y ■ q 2 ) 2 ) + N 2 (n , y)( qi • q 2 ) + N 2 (n , y)(y ■ qi )(y ■ q 2 ) + . . . (16) 



The expression ([14]) together with Eqs. (fT5|) and (fTBj) de- 
termine completely the gradient expansion of an arbi- 
trary correlation function F to second order in the den- 
sity gradients. These equations form the main result of 
the present work. If expansions to higher order gradients 
are required they can be derived without difficulty along 
the lines described above. To do this one needs to con- 
struct higher order symmetric polynomials in the wave 
vectors and carry out the required Fourier transforms. 
The consistency conditions Eq.© or equivalently Eq. ([T2)) 
then allow for a complete re-summation and leads to gra- 
dient coefficients depending on n(r). What is needed to 
determine the explicit form of the gradient coefficients in 
practice is the determination of the functions N m . How 
to do this for iV 1 and A^ 2 is described in the next section. 



C. Determination of N 1 and iV 2 

A practical calculation of the coefficients of the gradi- 
ent expansion of Eq. (|14[) requires the knowledge of the 
the functions N 1 and TV 2 . According to Eqs.Q and © 
these are defined as the first and second density deriva- 
tives of the correlation function F with respect to the 
density, i.e. 



<5F(r',r) 
Sn(r") 
<5F(r',r) 



5n{T")5n(r'"y no 

To derive useful expressions for these functions it is con- 
venient to transform derivatives with respect to the den- 



L = iVW-r) 

iV 2 (y,r"-r,r"'-r). 
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sity to derivatives with respect to the external potential 
as such derivatives appear naturally in perturbation the- 
ory. We can then use the well-developed tools of pertur- 
bation theory to calculate these functions. Let us there- 
fore define the new functions 



. 5 2 F(r',r) 



■F 2 (y,r"-r,r" 



6v(r")Sv{r"') 



,/,\ |n - 



(17) 



(18) 



Since we evaluate these functions for the homogeneous 
electron gas they only depend on differences between the 
spatial coordinates. It will therefore be convenient to also 
define their Fourier transforms by 



ip-ri+'iq-r 2 



t\( \ [ dpdq , 
(ri,r2) = J J2^ T (P ' q 

x-2/ \ f dkdpdq 2 

T (n.ra.rs) = J ^ 9 F (k,p,q) 

x ik-ri+ipr 2 +iq-r 3 

as well as their partial Fourier transforms 

^(y,q)=/^3^(p,q)e- y 

^ 2 (y,P,q)= / ^-F 2 (k,p,q)e^. 

The chain rule of differentiation gives a connection be- 
tween the derivatives with respect to the density and the 
derivatives with respect to the potential of the correlation 
function F. We have 

<5F( ri ,r 2 ) _ /^^MfaW 



Sv(r 3 ) 
<5 2 F(r 1; r 2 ) 
Sv(r 3 )Sv(r 4 ) 



dr 5 dr 6 



dr> 



<5n(r 4 ) Sv(r 3 ) 
JF(ri,r 2 ) 5 2 n(r 5 ) 
Sn(r 5 ) 5v(r 3 )Sv(r 4 ) 



<5 2 F(n,r 2 ) fa(r 5 )fa(r 6 ) 
8n(r 5 )Sn(r 6 ) Sv(r 3 ) Sv(r 4 )' 



(20) 



We see that in these expressions the first and second order 
derivative of the density with respect to the potential 
appear. We therefore define the linear density response 
function by \ 



<5n(ri) f dq 

x , L = X(r 2 - ri) = / — — jx(q 

ov{r 2 ) J (27r) d 



,»q-(.r2— ri ) 



as well as the second order density response function x 2 
by 



<5 2 ?i(ri) 
Sv{r 2 )Sv(r 3 ) 



lno = X 2 (r2 - ri,r 3 - n) 



[ rfprfq y 2 fr, a V*P-(r2-r 1 )+ t q.(r 3 -r 1 ) 

y (2tt)6 x ip ' qje 



Using these definitions we find by Fourier transforming 
Eqs.® and O that 



A^v.q) 



iV 2 (y,p,q) 



x(q) 



(21) 



[-^ 2 (y, p, q) - ^(y, p + q)x 2 (p, q)] 



x(p)x(q) 



(22) 



These equations give the desired relation between the 
density derivatives iV 1 and iV 2 and the potential deriva- 
tives J-" 1 and J-" 2 of the correlation function F. Relations 
between the higher order derivatives N m and T m can be 
derived in a completely analogous way. From Eqs. (|21[) 
and (f2"2"| we see that to calculate the functions TV 1 and 
./V 2 , and hence the gradient expansion coefficients of F 
to second order in the gradients, we need to calculate 
the density response and the change in F to second 
order in a perturbing potential 5v. The problem is 
thereby reduced to doing perturbation theory on the 
electron gas. This is a well-developed field of research. 
One option is to use diagrammatic perturbation theory 
[30| . In the following section we will use standard 
perturbation theory to obtain these response function 
for the case that F is the one-particle density matrix 
and use that to calculate the gradient expansion of the 
exchange hole. 



III. GRADIENT EXPANSION OF THE 
ONE-PARTICLE DENSITY MATRIX AND THE 
EXCHANGE HOLE 

A. The one-particle density matrix 

As an application of our formalism we carry out the 
gradient expansion of the one-particle density matrix for 
a noninteracting system with density n(r). This prob- 
lem has received large interest in the past since both the 
exchange energy E x [n] as well the Kohn-Sham kinetic 
energy T s [n] are explicit functionals of such a noninter- 
acting density matrix 0, Q . Therefore a gradient expan- 
sion of this density matrix directly leads to a gradient 
expansion of these two functionals. Such a gradient ex- 
pansion was first carried out by Gross and Dreizler [3l| 
using the Kirzhnits expansion 0, |32[ • This expansion is 
adjusted to the specific form of the noninteracting one- 
particle density matrix and can not be generalized to ar- 
bitrary correlation functions. The method presented in 
this work can, on the other hand, be applied to arbitrary 
correlation functions, but to demonstrate the soundness 
of our formalism we will use it to provide an alternative 
derivation of the result obtained earlier by Gross and 
Dreizler using the Kirzhnits expansion. An advantage of 
our derivation based on nonlinear response theory is that 
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it avoids the turning point problem encountered in the 
Kirzhnits approach [2|. 

Let us start by defining the one-particle density matrix 
in second quantization as 



7 



(r',r) = ^;(*|Vit( ra )^( r V)|*) 



where a is a spin coordinate. We will consider the case 
of spin-compensated systems. Then, in a noninteracting 
system with 27V electrons, the density matrix is given in 
terms of one-particle wave functions by 



N 

7 (r',r) = £2 W (r>*(r) 



(23) 



where the pre-factor 2 is a spin factor. The one-particle 
states are eigenstates of a one-particle Schrodingcr equa- 
tion 



1. 



V 2 +v(r) tpj(r) = ejtpjir) 



(24) 



where v(r) is the external potential. To determine the 
gradient expansion coefficients of 7 we need to determine 
how the density matrix changes if we make small changes 
v(r) — > v(r) + 5v(r) in the external potential. More pre- 
cisely, we need to calculate the functional derivatives 



7 1 (ri,r 2 ,r 3 ) 



7 2 (ri,r 2 ,r 3 ,r 4 ) 



fry(ri,r 2 ) 
5v(r 3 ) 
(5 2 7(r!,r 2 ) 
6v(r 3 )6v(r 4 ) 



which are the equivalents of the functions T 1 and T 2 of 
Eqs. (fl7| and (JTHJ) for the specific choice of correlation 
function F = 7 (note that they become functions of two 
and three independent vectors respectively when evalu- 
ated for a homogeneous system) . To do this it is sufficient 
to know the functional derivatives of the orbitals tfj and 
eigenvalues ej with respect to the potential v. These are 
simply given by doing first order perturbation theory on 
the one-particle Schrodingcr Eq. (|24)) . Wc find 



Sv(r) 

<Vj( r ) 
Sv(r') 



<^(rK(r')^(r' 



(25) 
(26) 



With these relations we can differentiate Eq. (f23"f twice 
with respect to the potential. Afterwards we then need 
to insert the plane wave states f-kiv) and their eigenval- 
ues ek appropriate for the homogeneous electron gas into 
the final expressions. In order not to interrupt the pre- 
sentation these calculations arc given in Appendix [Al and 
we simply present the results of these calculations here. 
If we define the Fermi factors by 

n k = 9(k F - |k|), 
where 9 is the Heaviside function and where k-p = 
(37r 2 no) 1 / 3 is the Fermi wave vector, then the resulting 
expressions are given by 



7 x (y>q) 



7 2 (y>p>q) 



where 



dp n 



p+q 



!E e ip-y 



(2tt) 3 £p +q -e p 
^-^e ik ^[$(k,k + p + q,k + p 



(27) 



$(k,k + p + q,k + q) 



$(k,p,q) = 



(28) 



(29) 



In this expression e k = |k| 2 /2 are the one-particle ener- 
gies. 

It remains to calculate the first and second order den- 
sity response functions %(q) and x 2 (p, q) to evaluate the 
functions N 1 and N 2 of Eqs. (JUJ) and P2]). Fortunately, 
for the specific choice F = 7 that we made these density 
response functions are simply given by 



X(q)=7 1 (y = 0,q) 

x 2 (p,q) = 7 x (y = o,p,q). 



(30) 
(31) 



They are therefore obtained directly from Eqs. (|27|) and 
(|2^|) so that we do not need to do any extra work to 
calculate them. 

To calculate the gradient coefficient to second order in 
the gradients we now need to expand 7 1 and 7 2 to second 
order in the wave vectors p and q. Since these calcula- 
tions are somewhat lengthy we do not present them here 
and refer to Appendix [B] for a description of the main 
steps. The resulting expansions are given by 
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7 1 (y,q) = 



7 2 (y,p,q) 



smz 



1 - 



my (q-y) : 



12fc| 



cos(z) 



i 1 

c °s(^) - ^(p + q) • ycos(z) + '- 



(P 2 + q 2 + P ■ q)( cos ( 2; ) + zsin(jz)) 



-^[(P • y) 2 + (q ■ y) 2 + 3((p + q) • y) 2 ] cos(z) 



(32) 



(33) 



where we defined z = fay. The expansions for the first and second order response functions x an d X 2 °f Eqs. (|30j) 
and (f31"j) are obtained by taking y = in these expressions. Together which Eqs. (|32|) and (|33|) these expression then 
immediately determine the functions TV 1 and N 2 from Eqs. (j2"Tj) and (f22|) . The final result of this calculation to second 
order in the wave vectors is given by 



^(y.q) 



^ 2 (y,P,q) 



smz 

z 



«qy 



q-y 



COs(z) 



2 6 

sin(z) 



q 2 sin(z) — z cos(z) 
12ftf z 



i- 2(p + q)-y - g[(p-y) 2 + (q-y) 2 : 



i 



+ TT77 2"(P +q +P-q)(3cos(z) + zsin(z)-3 



12/c 2 



sin(,z) . 1 . w N r sin(z) , . , 

—^ L ) + y^(p ■ y)(q • y)t 4 ^ - 3 



By comparison of these two expression to the expansions (|15|) and ()16|) we can directly read off the gradient coefficient 
functions of Eq. ([l"4"l) to be 



Nl(y) 



12fc£ 



[z 2 j (z)-3zj 1 (z)} Nl(y) 



12fcf 



[jo (2) + 3zji(z)] 



where now we replaced no — »■ n(r) and hence the Fermi wave vector fa = (37r 2 n(r)) 1 / 3 that appears in these expressions 
is a spatially varying function. We further introduced the spherical Besscl functions 



3o(z) 



sin z 



Ji(*) = 



sin z — z cos z 



By inserting the gradient coefficient functions into Eq. (jl4p we find the following explicit gradient expansion for the 
one-particle density matrix 



7(r',r) 



fc F JiO) 1 . 

— + ~3o(z)y ■ Vn r 

7T Z Z 



12*£ ZJ ' l(z)V n(r) + 6fcf 2 jo(z)( y ' V ) n ( r ) 



24fc| 



(« 2 io(«)-32iii(«))(Vn(r)) a 



24fc| 



(jo(z) + 3zji(z))(yVn(r)) 2 



(34) 



This expression is the main result of this section. After some manipulations we can rewrite it in an equivalent form 
as 

J_\ h3 3i(z) 1 V 2 fc 2 

r2 



7 (r',r) 



7T- 



Z 
2 \2 



24 fa 



1 1 

12 fc^ 



(v-(Vfc| -^ 2 j (z) + V 

V y / y 4 



1 (Vfc 2 ) 
96 fc| 



1 1 



32 fcp 



(j (z)z 2 - zj^z)) + -_^( Vfc 2 • ^ ) (z 2 j (z) - z A h{z)) 



(35) 



This expression becomes identical in form to that derived by Gross and Dreizler [31j using the Kirzhnits method 
after eliminating the effective Fermi vector of their work in favor of the density (this requires inversion of Eq.(19) of 
reference [3l| ). We close this section with a final remark on our result. We note that the gradient expansion to finite 
order breaks the symmetry 7(r,r') = 7(1"', r)*. However, the full gradient expansion as described by Eq. ([TTl) has this 
symmetry, which means that the symmetry is restored by taking all higher order gradients into account whenever the 
series converges. The symmetry violation has clearly consequences for the quantities that are derived from it, such as 
the exchange hole, as it introduces ambiguity in defining such functions. To be in accordance with existing literature 
we will in the next section adopt the definition of the exchange hole that was used by Perdew Q. 
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B. The exchange hole and energy 

We finally stress a few points on the properties of the exchange hole as derived from the gradient expansion and 
present an alternative derivation of the gradient expansion for the exchange energy compared to that of Dreizler and 
Gross which has the advantage of avoiding the regularization of divergent integrals [311 ] . In doing so we confirm a 
result obtained by Langreth and Perdew [ljj. The exchange hole can be calculated directly from the one-particle 
density matrix as 

| 7 (r + y, r)| 2 



Px(r + y|r) 



2n(r) 



Inserting the expansion (|34[) of for the one-particle density matrix and retaining terms to second order in the gradients 
yields the expression 

Px(r + y|r) = -|n(r)^ - |Mf^lM (y . Vn(r)) + _L Jl(z)2v 2 n(r) - i*jo(*)ii(*)£ . V ) 2 n(r) 

+ £ji{z)(zMz) ~ 3ji(z))(Vn(r)) 2 + ^(^Mfl + 3jl (z) 2 - 3j (z) 2 )(y • Vn(r)) 2 + . . . (36) 

F F 

As was pointed out by Perdew Q the second order gradient expansion of the exchange hole does not satisfy the 
negativity and sum rule constraints 

p x (r + y|r)<0 J dy /9 x (r + y|r) = -1. 

The violations of the negativity constraint is readily verified. However, the violation of the sum rule is not immediately 
obvious from Eq. (|36[) . The sum rule can, however, be conveniently analyzed in momentum space. To do this we define 
the Fourier transformed exchange hole to be 

Px (k|r)= J rfy Px (r + y|r) e -^ (37) 

This function has the form 

p x (k|r) =p°(k|r) + a x (k, n) k • Vn(r) + a 2 (k, n) V 2 n(r) + a 3 (k, n) (k • Vn(r)) 2 

+ a 4 (k, n)(Vn(r)) 2 + a 5 (k, n) (k ■ V) 2 n(r) + . . . (38) 

where we defined the unit vector k = k/fc and k = |k|. The explicit form of the functions p x and oij follow directly by 
Fourier transforming Eq. (|36[) and are given in Appendix [Cj The coefficient functions ctj are tempered distributions [33j 
which have mathematically well-defined Fourier transforms. As follows directly from Eq. (|37p the sum rule condition 
in momentum space translates to the requirement that p x (k = 0|r) = —1. Since p x (k = 0|r) = —1 (see Eq. (|Cl|) ) 
this implies that the sum rule would be fullfilled whenever aj(k = 0, n) = 0. However, as is clear from Eqs. (|C2[) - (|C6[) 
in the Appendix this is not satisfied for a\. whereas a.3 and a§ even diverge for k — > when interpreted locally as 
functions rather than distributions. The gradient expansion of the exchange hole does therefore indeed not satisfy the 
sum rule. These divergencies, however, cancel when we calculate the system averaged exchange hole in momentum 
space which is given by the expression 

N( Px (k)) = J drn(r)p x (k|r) 

r (p°(k)) +k- J dr?i(r)ai(£:,n(r))Vn(r) + ^ J drn(r)ftj(fc,n(r))ain(r)^n(r). (39) 



N( 



»J=1 ' 



where N is the number of electrons. In this expression term under the integral sign we defined the tensor 
(p x (k)) is the system average of p x (k|r) and in the last 



fcj(n, k) = a L (n, k)^f + a T (n, k){s tj 



k' 2 
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This tensor has a longitudinal part with coefficient 
and a transverse part with coefficient ax which describe 
the contributions to the system-averaged hole of density 
gradients Vn parallel and perpendicular to the momen- 
tum vector k. These coefficients are calculated from the 
functions a,- as 



«T =0:4 



1 d(na 2 ) 
n 



dn 
1 d(na 5 ) 
n dn 



it is now a simple matter to calculate the exchange energy 
from 

_ . . N f dk 4tt 

We insert Eq. (|3T)f and do the angular integrations first. 
It is therefore convenient to define the spherical and sys- 
tem averaged hole in momentum space by 

<0°*( fc )» = 4^ / d ^ {P*M) 



as a short calculation on the basis of Eq. ([38|) will show. 
From these equations and the explicit expressions for ctj 
given in Appendix [C] we can readily calculate the explicit 
expressions for ol,t- If we define the dimensionless vari- 
able k = fc/(2fcp) then we have 



such that 



ax 



'32nfc| 



32nfcf 



Z T (k) 



where the functions Zl,t have the explicit form 



Z T (x) = - 4x6(1 -x) + -S(x - 1) 



Z h ( x ) = - 4x9(1 -x) + S(x - 1) + -6'(x - 1) 

o 

These expressions are in accordance with the results of 
Langreth and Perdew (see Eq.(3.55) of reference [l9|) and 
this independent result therefore confirms the correctness 
of our alternative derivation. The interesting observation 
is now that 



lim Zl t(&) 

fe^O 







and hence /3y(fc,n) — > when k — > 0. This implies that 
the last term in Eq.(f3T)f does not contribute to the sum 
rule of the system averaged hole. The same conclusion 
can be derived for the second term in Eq. (f3T)]) . Since 
ot\(k = 0,n) = i/(4nkp) (see Eq. (|C2[) ) is a local func- 
tion of the density the integrand of the first integral in 
Eq. (|3l?)) is a total derivative and vanishes (assuming ei- 
ther vanishing densities at infinity or periodic boundary 
conditions). We therefore find that 

lim( Px (k)) = hm(p°(k)) = -1. 

fc— s-0 k— >0 

This implies that system averaged exchange hole ob- 
tained from the second order gradient expansion does 
fulfill the sum rule, although the exchange hole itself docs 
not. We note, however, that this is only achieved by in- 
tegrating over both positive and negative contributions. 
When the positive contributions to the exchange hole are 
removed the sum rule is only recovered for a finite hole ra- 
dius [j| [l(| ■ We finally note that with expression Eq. ([55)1 



E K [n] = — f dk(( Px (k))). 

7T JO 



(40) 



From Eq.© we then find that 

N(( Px (k))) = N((p°(k))) 

f 1 2 

+ / dxn(r)[-a L (k,n) + -a T (k,n)](\7n) 2 

= N{{p° x (k))) + J dr 3^*(*)(Vn) 2 , (41) 



where we defined 



Z x (k) = ^ L (fc) + \z T (k). 



(42) 



The exchange energy is obtained by inserting the expres- 
sion for the averaged hole of Eq. (|4*T|) into Eq.(T5D]). This 
then finally gives the expression 



E x [n] = - ^(^) 3 J drn3(r) 

+ / drB x (n(r))(Vn(r)) 2 



(43) 



where 



B x (n) 



164 JO 



dk Z x (k) 



7 



4327r(37r 2 ) 1 /3 n 4/ 3 



and where to calculate the first term in Eq. (|4"3"]) we used 
Eq. ([Cl[) . The coefficient i? x is the same gradient coeffi- 



cient as obtained by Gross and Dreizler [31| and earlier 
by Sham [34[. However, the correct analytic exchange 
gradient coefficient is known (35l - [39j to be a factor 10/7 
larger. The reason for this discrepancy is clearly de- 
scribed by Svendsen and von Barth [38[ who showed that 
the Coulomb interaction is too singular to allow for the 
interchange of the operations of integration and the ex- 
pansion in wave vectors. The problem does, for instance, 
not appear when Yukawa-screened Coulomb interactions 
are used [38|. The conclusion is that there is nothing 
wrong with the Kirzhnits method, or the nonlinear re- 
sponse theory derivation of the gradient expansion of the 
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exchange hole presented here, but one has be aware that 
for Coulomb interactions the procedures of expanding 
correlation functions in terms of density gradients fol- 
lowed by integrations involving Coulomb potentials may 
not yield the same result as directly expanding the inte- 
grated quantity in terms of density gradients. For this 
reason the original GGA of Pcrdew and Wang [l(| based 
on the gradient expansion of the exchange hole was later 
reparametrized [lj| to accomodate the correct gradi- 
ent coefficient for the exchange energy. 

IV. CONCLUSIONS AND OUTLOOK 

We derived a general scheme based on nonlinear re- 
sponse theory to calculate the density gradient expan- 
sion of general correlation functions and showed that in 
order to express the gradient coefficients in terms of the 
full density profile summations to infinite order must be 
carried out over response functions of arbitrarily large or- 
der. A consistency condition was derived to do this. The 
formalism was used to derive the gradient expansion of 
the one-particle density matrix and the exchange hole to 
second order in the gradients. We confirm the derivation 
of Dreizler and Gross that used the Kirzhnits expansion 
method. We further analyzed the exchange hole in mo- 
mentum space to derive that the system averaged hole to 
second order in the gradients satisfies the sum rule and 
to derive the gradient expansion of the exchange energy 
without the need to regularize divergent integrals. 

The scheme that we presented is very general and can 
be applied to more general correlation functions. A nat- 
ural application of the scheme would be to calculate the 
gradient expansion of the correlation hole p c . Regarding 
the gradient terms of the correlation hole little is known. 
We essentially only have some detailed information on 
the long-range properties of the spherically and system 
averaged hole. This information comes from the work of 
Langreth and Perdcw who calculated ((p c (k))) within the 
RPA. In our notation their results (see also Appendix C 
of [IH) can be summarized as 

N(( Pc (k))) = N{{p° c (k))) + J dv^jz c {n,k){Vnf 

where the function z c has been parametrized by Langreth 
and Mehl @, HJ as 

Zc{n , k ) = ^l e -^£ (44) 

where fc s = (Ak-p/ir) 1 ^ 2 is the Thomas-Fermi or screening 
wave vector. We can transform to real space to obtain 
the following expression for the spherically and system 
averaged correlation hole 

W<Ms/)» = N((p°(y))) 

+ /*<i^(TT(J)W <V " )2 <45) 



We see from Eq. (|44)) that z c (n, k = 0) ^ and as a con- 
sequence the sum rule ((/f c (0))) = for the correlation 
hole is not satisfied. We know, however, that the RPA 
becomes accurate in the high density limit and since from 
Eq. (|45|) the gradient coefficient of the correlation hole de- 
pends on k s y we see that high densities are equivalent to 
large separations y. Similarly, low density corresponds 
to short-distance behavior. However, RPA can not be 
trusted in this region and we have no precise knowledge 
on the small distance behavior of the gradient coefficient 
of ((p c (y))). A model for the general short-range behav- 
ior was proposed [l6hl8j on physical arguments (it should 
not affect the Coulomb cusp of and the on-top value of 
the LDA hole) after which the real-space cutoff proce- 
dure was applied (see e.g. Fig. 4 in [17| and Fig. 5 in 
[lj| ) to obtain a GGA for correlation. It would, how- 
ever, be desirable to have a first principles approach to 
the calculation of the short range properties of the gra- 
dient coefficient of the correlation hole. As follows from 
our derivation this requires the knowledge of the func- 
tions (|17p and (JTSJ) , or at least their expansion to second 
order in wave vectors, for the case that F represents the 
pair correlation function or the correlation hole of the 
electron gas. It may therefore be worthwhile to use our 
current scheme to explore these response functions be- 
yond the RPA. For short range correlation an approach 
based on diagrammatic summation of ladder diagrams 
suggests itself. Of course, for the development of density 
functionals for general systems beyond the weakly inho- 
mogencous regime it is not sufficient to use the straight- 
forward gradient expansion [23j . However, general short- 
range features, such as the way the exchange-correlation 
hole deforms close to the reference electron in the pres- 
ence of a density gradient can be transferred to more 
general systems than the weakly inhomogeneous ones. 
Perhaps in combination with truly nonlocal functionals 
for the exchange-correlation hole [4LJ this can lead to the 
development of more accurate density functionals. This 
will be a topic of our future research. 

Finally, we would like to make some general remarks 
on the extension of our derivations to temperature- or 
time-dependent systems. In the case of temperature- 
dependent systems the expectation values of observables 
are traces over a grand canonical ensemble. By Mermin's 
theorem [42[ such expectation values are still functionals 
of the density and therefore the derivations in this work 
can be carried out unchanged. In fact, the dependence 
on temperature is probably going to improve the conver- 
gence properties of the gradient series by smoothening 
of the sharp Fermi functions which, for instance, were 
the cause of the oscillatory nature of the gradient coeffi- 
cients of the one-particle density matrix. The situation 
of time-dependent systems is more delicate. Although 
the existence of a density-potential mapping has been 
well-established [43l - |4r| and hence a density functional 
theory can be formulated, the time-variable induces se- 
vere complications. As was shown by Vignale and Kohn 
[47l [HI the temporal and spatial non-locality of time- 
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dependent density functionals are intimately correlated. 
Any temporal non-locality (or frequency dependence in 
the language of equilibrium response functions) induces 
long-range spatial properties that prevents the gradient 
expansion from existing. A way out of this problem is 
to use new variables such as the current-density [53, EH 
or the local density deformation density in a Lagrangian 
frame (49l - [5l| for which a local density approximation 
and a corresponding gradient expansion does exist. This 
has been exploited in Ref.[52j to construct a GGA func- 
tional within time-dependent current density functional 
theory. The study of correlation functions in the same 
fashion remains a challenge for the future. 
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Appendix A: Calculation of 7 1 and r y 2 

In this Appendix we will derive the expression for the 
first and second derivative of the one-particle density ma- 
trix 7 with respect to the potential v. For a clearer inter- 
pretation and compactness of notation we will only put 
a comma between the variables representing the original 
coordinates of the density matrix and the variables that 
appear as argument of the potential variations. Direct 
differentiation of Eq. ([2"B")) using Eq. (|2l))) then gives 



7 1 (23,1) = 



M23) 
Sv(l) 



4^ <P,-(1W?(1W2)<^(3) 

2^-/0 . ( A1 ) 



where we used the short notation j = rj and introduced 
the occupation numbers fj = 1 for an occupied state and 
fj = for an unoccupied state. Inserting plane wave 
states appropriate for the electron gas we find that 



7 1 (r 2 r 3 ,n) = 2 



dqdp re p+q 
(2^) 6 e P+q 



P gip'(r 2 — r 3 )-Mq-(ri— r 3 ) 



We therefore find that the function 7 1 has the simple 
form 



7 



1 (p,q) = 2 



'■p+q 



tp+q 



in Fourier space. This is precisely the integrand of 
Eq.(f2~T|). The function 7 2 can be calculated by straight- 
forward differentiation of Eq. (|Aip with respect to the po- 
tential 



7 2 (23,14) 



<V(23, 1) 
(5w(4) ' 



To obtain the explicit form of this function we have to dif- 
ferentiate both the orbitals and the eigenvalues using Eqs. 
([251) and ([26|) . Let us denote the term that arises from 
the differentiating the orbitals by X and the term that 
arises from differentiating the eigenvalues by Y. Then we 
find that 



7 2 (23, 14) = A(23, 14) + Y(23, 14), 



where 



X(23,14) = 2j>i(2)^(3) 

i,j,k 

4>(z^V fc (lM(lM(4)^(4) 
+ $(jik)<p j (lM(l)<p*(4) V > k (4) 

OO 

y(23,i4) = -25^(2)^(3) 

i,3 



(A2) 



$(iii)(| W (4)| a -| W (4)| 3 )^(lM(l) (A3) 



and where we further defined 



fk fi fj fi 



(A4) 



The function $ has the useful properties = and 

$(ijfc) = The function 7 2 (23, 14) = 7 2 (23,41) is 

symmetric in the indices 4 and 1 due to the fact that the 
differentiations with respect to the potential commute. 
This is, however, not obvious from Eqs.(jA2| and (|A3[) . 
We therefore want to rewrite the form of the function 7 2 
in such a way that this symmetry is obvious. To do this 
we first note that since $>(iji) = the terms with 

k = i and k = j in Eq. ([A2[) sum up to a function Z of a 
similar form as Y in Eq. (|A3p . namely 



Z(23, 14) = -2 X> (2)^(3) 

i-.j 

$^)(|^(1)| 2 -|^(1)| 2 V,(4)^(4). (A5) 
We can therefore write 7 2 (23, 14) as 

oo 

7 2 (23,14) = 2 Vi(2)^(3) 



$(jzfc)^(l)^(l)^(4)^(4) 
y (23, 14) + ^(23, 14) 



(A6) 



It is easily seen that the sum of Y and Z is symmetric 
under interchange of 1 and 4. However, this is not obvi- 
ous in the first term of the equation above since at first 
sight Q(jik) does not appear to be equal to §(ijk) since 



<f>{jik) 



1 



fk fj fi fj 



(A7) 
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which seems to be different from expression (|A4[) . How- 
ever, this is just appearance. The reason is that the oc- 
cupations can only attain the values zero and one. For 
the case /, = fj = 1 or /, = fj = we see directly that 
Eq. ([A4[) and (|A7|) attain the same value. A little inspec- 
tion shows that this is also true for cases fi = 0, fj = 1 
and /, = l,fj = 0. We therefore find for k ^ i,j that 
$>(ijk) = $(jifc). Therefore Eg. (|A6|) can be simplified to 

oo 

7 2 (23,14) = 2 Vi( 2 M( 3 ) 

ijk,k^{i,j) 

<S>(ijk) [^(1)^(1)^(4)^(4) 

+ ^•(1)^(1)^(4)^(4)" 
+ r(23,14) + Z(23,14) 



(A8) 



This expression is now explicitly symmetric in the indices 
1 and 4. Let us now evaluate this expression for the 
homogeneous electron gas. In the electron gas the one- 
particle states are plane waves 

Vk(r) = 



V 



where V is the volume of the system. Since |</>k| 2 = 1/V 
it follows from Eqs. (|A3| ) and (|A5|) that the terms Y and Z 
are identically zero and the function 7(23, 14) of Eq. (|A8[) 
therefore attains the form 



7 2 (23,14) 



2 

y 3 



$(k,p,q)e l(kr2 - pr3) 

k,p,q,q#(k,p) 



gi(q-k)-ri+i(p-q)-r 4 , g j(p-q)-ri+i(q-k)-r 4 



where we defined 
$(k,p,q) = 



rig - n k _ Up - n k 
£ q — £k e P — £k 



In this expression n p = 0(e F — e p ) is the zero-temperature 
Fermi function and £f = k F /2 is the Fermi energy. When 
we replace the sum by an integral the restriction q^k,p 
is a region of measure zero and we can therefore freely 
integrate over all wave vectors. After a few substitutions 
we obtain 



/ dkdpdq 2 

1>F ' 



7 2 (23,14) 



7 2 (k,p,q) 



gik-(r2— r3)+ip-(n— r3)-Hq-(r4— r3) 

where we defined 

7 2 (k,p,q) = 2[$(k,k+p+q,k+p)+$(k,k+p+q,k+q)] 
This gives the integrand of Eq. 



Appendix B: Expansion of 7 1 and j 2 

1. Expansion of 7 1 

To determine 7 1 to second order in the wave vectors 
we have to expand the function 



^•<»= 2 /(! 



dp n 



)3 e 



p+q n p gip-y 
P+q _ e P 



(Bl) 



to second order in powers of q. This can be done by 
expanding the integrand in powers of q. If we denote 
A = e p + q — e p = p q + <7 2 /2 with q = |q|, then we 
can expand the Fermi function n P + q in a distributional 
Taylor series as 



'p+q 



dn A 2 d 2 n A 3 d 3 n 



Since n p = 8(e F — e p ) is the Heaviside function we have 

A 2 

"p+q = rip - AS(e F - e p ) + — 5 (e F - e p ) 



A 3 



5"(e F - e p ) 



Inserting this into Eq. (|Bll) then gives 
7l(y ' q) = - 2 /(2^^ F - £p)elP ' y 



(B2) 



We see that in these integrals several derivatives of the 
delta function appear. We now use the general mathe- 
matical expression 



1 d 
ly~dx 

dx , 



E 



5(x - Xj) 



where the sum runs over all zeros of the function y(x), i.e. 
y(xi) = 0. Using this equation in spherical coordinates 
with y(p) = (fcp — p 2 )/2 we find 



S(e F 


-e P ) 


S'(e F 


-e P ) 


S"(e F 


-e P ) 


S"'(e F 


-e P ) 



1 

k F 



-Ls'( P -k F ) 

pk F 

5 (p- fcp) 8 (p- k F ) 



k F p 3 



k F p 2 



(B3) 
(B4) 

(B5) 



1 

k F 



1 /// 3 // 

-r$ (p-k F ) jS (p-k F ) 

p.1 pi 



—5 (p - k F ) 



(B6) 
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where in Eq. (|B6[) we also evaluated the third derivative of 
the delta function as we will need it later in the expansion 
of 7 2 . With Eqs. (|B3|) - (|B5j) we can readily evaluate the 
three integrals in Eq. (|B2[) . Let us denote these integrals 
by A, B and C. Then we find for the first integral 



A = - — 



2 

fe F 



dp 



5(p — k F ) e 



p y 



fcp sin z 



where z = k F y. The second integral gives 

B = J ^3(p-q + g 2 /2)5'(e F -e p )e^ 
= (q 2 /2 - • V y ) J j^8'(e F - Cp ) e*'* 



(q 2 /2 - iq • V y 

,,2 



1 



47r 2 fc T 



cos z + 



27T 2 fc F 

iq ■ y sin z 



cos(k F y) 



2tt 5 



-fe F - 



where in the second step we used Eq. (|B4[) . It therefore 
remains to calculate the last integral 



C 



dp 



(p-q+ g 2 /2)V(e F -ep)e'>" 



(2tt) 3 

However, up to order q 2 it is suffices to calculate 

c 1 r dp 



3./ (2tt) s 



-(p. q fS"(e F -e p )e^ 



3 V " "J (2tt) 3 
"g(q- V y ) 2 ^3-(cos(z) + zsin(0)) 



<7 2 , n , (q-y) 2 , sin(z) 

cos(z) H „ k F - 



67T 2 fc F 



6tt 2 



where in the second step we used Eg. (|B5|) . Adding the 
results 7 1 = A + B + C gives the expression (|32l) . 



We therefore only need to expand A(y, p, q) in powers 
of the wave vectors and symmetrize with respect to p 
and q afterwards. To do this we first need to expand the 
function 



$(k, k + p + q, k + p) 



£k+p — £k+p+q 

If we denote 



»k+p - "k _ "-k+p+q - »k 
£k+p — £k Ck+p+q — £k 



Ai = e k + P - £k = y + P • k 

a (p + q) 2 / \ , 

A 2 = e k+p+q - e k = h (p + q) • k, 



we find by expanding the Fermi functions in a distribu- 
tional Taylor series that 



$(k,k + p + q,k + p) = -S'(e F -e k ) 

--(Ai+A 2 )5"(e F -e k ) 
6 

+ —(A 2 + AxA 2 + A 2 )<f (e F - e k ) + 



With this expression we find that the function A(y, p, q) 
can be calculated as 



^(y,p,q) 

dk 



5 (e F - e k ) e 



iky 



(2tt) 3 

V(2^' 5 " (eF ~ ek)(Al+A2)elk ' y 



3 7 (2tt) 

i r dk 



12 J (2tt) : 



S"'(e F - e k )(A 2 + AiA 2 + A 2 ) e* ky 



2. Expansion of 7 

Analogously to 7 1 we can expand 7 2 (y, p, q) of Eq.(| 
in powers of p and q. To reduce our effort we note that 
we can write 



7 2 (y, p, q) = My, p. q) + A (y, q, p) 



where 



The three integrals appearing on the left hand side of this 
expression are sufficient to extract the powers to second 
order in p and q. Let us call these three integrals A\, A2 
and A3 in order of appearance. The first integral gives 



W(2^' (eF - £k)elky = 



2n 2 k F ' 



/dk 
- — — e lk y $(k, k + p + q, k + p) The second integral is given by 

(2-kY 



16 



A, 



i r dk 



3 7 (2tt) 3 

p 2 + (p + q) 2 

6 

p 2 + (p + q) 2 

6 



<5 (ep - £k) 



+ (p + q) 2 



(2p + q) • k 



,-(2p + q)-V y 



,-(2p + q)-V y 



cos(z) + z sin(z) 



-27T 2 kl 



cos(z) + zsin(z) 2 2 i 

= V 2 , 2fc| 1 V + (p + q) 2 ) - ^(2p + q) • y «»(*) 

where we used Eg. (|B5p in the second step. Finally the third integral has the explicit form 



-4, 



1 



dk 



+ ( 



12 J (2tt) 3 

(p + q) 2 , 



.p 



v 



• (P + q) 2 



(£- + p.k) 2 + (^-+p.k)( 
(p + q)-k) 2 lr(e F - ek )e lk - y 



+ (p + q) • k) 



However, we only need the terms to second order in p and q. Therefore, it is sufficient to calculate 



dk 



12 J (2tt) 3 



[(p ■ k) 2 + (p • k)((p + q) ■ k) + ((p + q) ■ k) 2 ] 6 (cf ~ eJ k y 

f dk a, 

12 [(p ' y y) 2 + (p ' y y)((P + q) ' y y) + « p + 1) ' v y) 2 ] / 7^3 5 e * ky 



= t( p ' Vy)2 + (P ' Vy)((p + q) ' Vy) + ((P + q) ' Vy)2 ] 



(27T) 

3 cos(z) — z 2 cos(z) + 3z sin(z) 



27T 2 fcf 

where to evaluate the integral over the delta function we used Eq. (|B6j) . To perform the derivatives in the last term 
we use that for an arbitrary function /(z) 



( p .V y )(q-V y )/(z) = (p.y)(q.y)^| 



d 2 / _ 1# 

dz 2 z dz 



+ (p ■ q) 



k¥_df_ 

y dz 



With this expression we find that 



As = ~ 247r 2 fc 3 (P 2 + P • (P + q) + ^P + q) 2 )( c os(^) + z sin (^)) 



1 



-[(p • y) 2 + (p • y)((p + q) • y) + ((p + q) • y) 2 ] cos(z) 



247T 2 fc F 

Collecting our results A = A% + A 2 + A3 we therefore obtain the following expression for 7 2 (y, p, q) 
7 2 (y,P,q) = 4(y,p,q) +A(y,q,p) 



1 



7r 2 fct 



cos(z) - -(p + q) • ycos(z) 



1 



12fc 2 



(p 2 + q 2 + p ■ q)(cos(z) + zsin(z)) 



-^[(P ' y) 2 + (q • y) 2 + 3((p + q) • y) 2 ] cos(z) 



This is equal to expression ([33 



Appendix C: Exchange hole in momentum space 

Here we present the explicit expressions for the coefficient functions a.j of Eq. (|38[) which are calculated by Fourier 
transforming Eg. (1361) . If we define the dimcnsionless quantity k = |k|/(2fcp(r)) then p° is given by 



-i + \k~ l -¥)e{\-k) 



(Cl) 
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which is simply the Fourier transform of the LDA exchange hole. We see that p x (0|r) = —1 and therefore the LDA 
hole satisfies the sum rule. The functions etj are given by 

ai =4^(1 - *) (C2) 



k 



+ + (C4) 



72n n,- 



"^wr^ 1 ^^^ -1 ^ (C5) 
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